Improved effective-one-body description of coalescing nonspinning black-hole binaries 
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We improve the efTective-one-body (EOB) description of nonspinning coalescing black hole bina- 
ries by incorporating several recent analytical advances, notably: (i) logarithmic contributions to 
the conservative dynamics; (ii) resummed horizon-absorption contribution to the orbital angular 
momentum loss; and (iii) a specific radial component of the radiation reaction force implied by 
consistency with the azimuthal one. We then complete this analytically improved EOB model by 
comparing it to accurate numerical relativity (NR) simulations performed by the Caltech-Cornell- 
CITA group for mass ratios q = (1, 2, 3, 4, 6). In particular, the comparison to NR data allows us to 
determine with high-accuracy (~ 10"'') the value of the main EOB radial potential: A{u; v), where 
u — GM / {R(?) is the inter-body gravitational potential and u = q/{q + 1)^ is the symmetric mass 
ratio. We introduce a new technique for extracting from NR data an intrinsic measure of the phase 
evolution, {Qui{ijj) diagnostics). Aligning the NR-completed EOB quadrupolar waveform and the 
NR one at low frequencies, we find that they keep agreeing (in phase and amplitude) within the NR 
uncertainties throughout the evolution for all mass ratios considered. We also find good agreement 
for several subdominant multipoles without having to introduce and tune any extra parameters. 
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I. INTRODUCTION 

The Effective One Body (EOB) formalism [1-5] has 
been proposed as a new analytical method for describing 
the motion and radiation of coalescing black hole bina- 
ries. One of its main aims is to provide analytical ^ grav- 
itational waves (GW) templates covering the full coales- 
cence process, from early inspiral to ringdown, passing 
through late inspiral, plunge and merger. The definition 
of the EOB formalism mainly relics on two sources of 
information: 

(i) high-order results of post-Newtonian (PN) theory; 

(ii) high-accuracy results from Numerical Relativ- 
ity (NR) simulations of coalescing black hole binaries 
(both in the comparable- mass case, ly = C(l), and 
in the extreme-mass-ratio limit, v 1). [Here, v = 
toi7712/(toi -I- 7712)^ denotes the symmetric mass ratio.] 

In addition, EOB theory has recently tapped useful 
information out of Gravitational Self Force (GSF) com- 
putations at order 0{i'). All this information is not used 
in its original form, but rather as a way to determine, 
or at least constrain, the structure of the few basic func- 
tions that enter the definition of the EOB formalism. For 
recent general reviews of the EOB formalism and its his- 
torical roots, see [6, 7]. 

The EOB formalism has been developed in a sequence 
of papers, both for nonspinning black hole binaries [1- 
3, 5, 8, 9] and for spinning ones [4, 10-14]. In addition, 
it has been extended to the case of tidally interacting 



neutron star binaries [15, 16]. For all those types of sys- 
tems, many comparisons between the predictions of EOB 
theory and the results of NR simulations have been per- 
formed [9, 17-29] and have demonstrated that it is pos- 
sible to devise accurate EOB waveforms by combining 
improved resummation methods [5, 8, 9], high-order PN 
results (see [30] for a review), and some nonperturba- 
tive information coming from high-accuracy NR results. 
These EOB waveforms can be used both in GW detec- 
tion and in GW parameter-estimation protocols. The 
EOB formalism can thereby crucially help detecting the 
GW's emitted by coalescing black hole binaries, since 
many thousands of waveform templates need to be com- 
puted to extract the signal from the broad-band noise, an 
impossible task for NR alone. The EOB formalism might 
also be crucial in allowing one to extract information on 
the equation of state of nuclear matter from observations 
of coalescing neutron star binaries [31]. An early ver- 
sion of the EOB waveform [28] has already been incor- 
porated^, and used [32] in the LIGO and Virgo search 
pipeline. 

In addition, some recent comparisons between NR 
studies of the dynamics of black hole binaries and its EOB 
description, have directly confirmed the ability of EOB 
theory to accurately describe several (gauge-invariant) 
aspects of the conservative dynamics of binary systems, 
such as periastron precession [33] and the relation be- 
tween energy and angular momentum [26]. 

The aim of the present paper is to improve the defini- 
tion of some of the basic elements of the EOB formalism 
both by including for the first time recently obtained an- 
alytical information, and by extracting, in a new way. 



^ Here we use the adjective "analytical" (instead of "semi- 
analytical") for methods that are based on solving analytically 



given ordinary differential equations, 
numerical tools to solve them. 



even if one needs to use 



See https:/ /www. lsc-group.phys.uwm.edu/daswg/projccts/lalsuite. html. 
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nonperturbative information from accurate NR simula- 
tions performed by the Caltech-Cornell-CITA group [34]. 
Though our study will be limited to nonspinning binaries, 
the EOB structures we shall improve (such as the basic 
EOB radial potential A{R)) are central, and must then 
be included both in the spinning and tidal extensions of 
the EOB formalism. 

The recent analytical progresses that we shall incorpo- 
rate here in EOB theory are: 

(i) 4PN and 5PN logarithmic contributions to the con- 
servative dynamics [35-38]; 

(ii) the 0{i') 4PN nonlogarithmic contribution to the 
conservative dynamics [36,38-40]; 

(iii) resummed horizon-absorption contributions to an- 
gular momentum loss [41, 42]; 

(iv) the radial component of the radiation reaction force 
implied by consistency with the azimuthal one [43] ; 

(v) an additional 3.5PN contribution to the phase of 
the (factorized [5, 8, 9]) quadrupolar waveform [44]. 

In addition, we shall bring up some novelties in the def- 
inition of the EOB formalism, and in the way to extract 
information from (comparable-mass) NR data. Namely: 

(a) we introduce a Pade resummation of the additional 
tail phases Sim of the factorized EOB waveform; 

(b) we show how to accurately extract from NR data 
the Qa;(w) function measuring, in a intrinsic way, 
the phase evolution of the (curvature) quadrupolar 
waveform; 

(c) we introduce a new way to improve the EOB wave- 
form during plunge and merger by matching it to 
the NR one at a specifically chosen (:^-dependent) 
NR time t^^^{i') around merger. More precisely, we 
impose [by using six next-to-quasi-circular (NQC) 
parameters] a contact condition between the 
amplitudes and the frequencies of the NR and EOB 
waveforms at an NR instant t^^j.{i'), correspond- 
ing to the maximum of the EOB orbital frequency 

iEOB 
'■Hpcak- 

The paper is organized as follows. In Sec. II we present, 
in a self contained manner, the detailed definition of 
our improved EOB formalism. Section III explains 
how to extract the (3a;(a;) function from NR data while 
Sec. IV revisits the extreme-mass-ratio case. In Sec. V 
we then complete our new EOB formalism by compari- 
son with several comparable-mass simulations performed 
by the Caltech- Cornell- CITA group. Section VI studies 
the structure of the main EOB radial potential {A{u) 
function) obtained from the latter NR comparison and 
Sec. VII discusses how to compute EOB waveforms for 
arbitrary values of v. We summarize our main conclu- 
sions in Sec. VIII, while some supplementary material 
is presented in several Appendices. In particular. Ap- 
pendix D gives the explicit expressions of the pim and 
dim bricks of the EOB factorized waveform we use. 



II. EFFECTIVE-ONE-BODY ANALYTICAL 
FRAMEWORK 

In this Section we shall present in detail the definition 
of the new (non spinning) EOB formalism, incorporating 
several recent analytical improvements that we shall use 
in this paper. Our presentation will be self-contained so 
as to allow readers to generate for themselves all our EOB 
results. We also intend to make available soon a public 
version of our EOB codes. 

The EOB formalism is made of three basic building 
blocks: (i) a EOB Hamiltonian that rcsums the conserva- 
tive two-body dynamics; (ii) a resummed EOB radiation 
reaction force that completes the conservative dynamics 
by causing the system to inspiral down to merger, and 
(iii) a resummed EOB inspiral-plus-plunge waveform, to- 
gether with a prescription for extending the waveform 
through merger and ringdown. Each one of these build- 
ing blocks has been developed in previous papers. In 
particular, the construction of the EOB Hamiltonian was 
initiated in Refs. [1, 3], while the definition of the re- 
summed, factorized inspiral waveform was initiated in 
Refs. [5, 8, 9]. Here we bring new (recently derived) the- 
oretical improvements to each element of the formalism, 
namely: (i) we include logarithmic contributions [35-38] 
to the EOB Hamiltonian; (ii) we include the effect of a 
resummed version of horizon absorption [41, 42] in the 
radiation reaction; (iii) we add a recently derived [43] 
radial component of radiation reaction; (iv) we include 
the 3.5PN contribution [44] to the phase 622 of the fac- 
torized quadrupolar waveform; (v) we resum S22, as well 
as some higher- multipoles Sgm^s, by Pade methods. All 
these improvements either add some new physics that 
was not included in the previous EOB models [22, 28], 
or improve [in the case of (v)] the robustness of the EOB 
resummations. We shall discuss them in detail in the 
sections below. 



A. Improved Hamiltonian: logarithmic 
contributions to the A function 

The conservative (non spinning) two-body dynamics 
is described, within the EOB formalism, by a Hamilto- 
nian H-eob{Q^ T Pi) , describing the relative motion Q' = 
Qi ~ Qh of the binary, and depending on two radial func- 
tions, A{R) and B{R), where R = \Q'^\ is the binary 
separation (in EOB coordinates). We are using phase 
space variables (i?, Pr, </?, P,^) associated to polar coor- 
dinates in the equatorial plane 9 = tt/2. Actually it is 
useful to replace the radial momentum by the mo- 
mentum Pfj^ = {A/ B)^/"^ Pr conjugate to the "tortoise" 
radial coordinate R^ = J dR{B / A)^/"^ . Furthermore, it 
is convenient to use suitably rescaled dimensionless vari- 
ables: 

R _ Pr, Py T 

"■"GM'^"*" u'P^'aGM' GM- ^> 
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Here, and in the following, we use the notation 

,,r , _ mmi2 _ _m2 

M ^ mi + 7712, fJ- = ■ , ^= T7^ 1= ■ (2) 

?Tll + 7712 M TOl 

Note that the dimensionless symmetric mass ratio v = 
mim2/(TOi+m2)^ = 9/(9+1)^ varies between (extreme 
mass-ratio case) and j (equal- mass case), and that we 
shall conventionally consider that toi < m2, so that q > 
1. In addition we generally set c = 1, and shall also often 
set G = 1 in the following. 

With the above notation, the /z-rescaled (real) EOB 
Hamiltonian reads 

HEOB{r,Pr,,Pip) = = + 21^ {Hcff ~ 1) , (3) 

where H^a denotes the (/i-rescaled) effective EOB Hamil- 
tonian, given by 



^(r)( 1-1-^ + ^3^ 



(4) 



with Z3 = (4 — Zv). 

The (rescaled) EOB Hamiltonian (3) leads to equations 
of motion for {r,ip,Pr,,Pip) with respect to the rescaled 
time t = T/GM, Eq. (1) , of the form 



^ = O - ^-^EOB 
dt dp^ 

dr /A\^^^ dH, 



'EOB 



dt \B J dpr, 
dp^ 



dt 



dt \B 
which explicitly read 

dt ~ 



2EOB 



dr 



(5a) 

(5b) 
(5c) 
(5d) 




where A! = dA/dr. In these equations, = T/n denotes 
the ^-rescaled radiation-reaction force. Its explicit form 
will be given in Sec. II D below. 

Let us now define the explicit forms of the two basic 
EOB radial functions A{r) and B{r) entering the Hamil- 
tonian (3). One of the main theoretical novelties of the 



EOB model used in the present work is the inclusion in 
A{r) (which plays the role of the main radial potential 
in the EOB Hamiltonian) of the recently computed log- 
arithmic contributions appearing at the 4PN and 5PN 
levels [35-38]. If we first focus on the Taylor- expanded 
version of the A potential, it has, when considered at the 
5PN level, the form 

^Taylor ^ I ^ 2u + 2vU^ + (— - VU^ 

y 3 32 j 
V \_a%{y) + CL^{v) Inu] + v \a%{v) -V a^l\v) Inw] ^t^ 

(7) 

where u = GM/R =l/r denotes the (EOB) dimension- 
less gravitational potential, and where 



64 
7004 



144 



105 



(8) 
(9) 



denote the analytically known logarithmic contributions, 
while ag(i') and a,Q{v) represent currently uknown, non 
logarithmic j/-dcpendent 4PN and 5PN contributions to 
A{u). Following the EOB methodology initiated in Ref. 
[3], we do not use the Taylor-expanded radial potential 
^Taylor to define the EOB Hamiltonian, but use in- 
stead its (1,5) Pade approximant, namely 



A{u- al{v),al{u)- u) ^ P}AA^^y'°\u)] 
1 -I- niu 



1 + + i 



(10) 



where the coefficients ni and d.j appearing in the numera- 
tor and the denominator of the Pade approximant depend 
rationally on 05,0^,1^ and Inu. As is well known, Pade 
approximants can sometimes exhibit "spurious poles" in 
u. In the case we shall investigate here (with a§ fixed to 
the value in Eq. (11)) such a spurious pole starts appear- 
ing when Og < —131. This will not affect our analysis 
because we shall work in the range > —123. 

The logarithmic-dependent 5PN-Pade-rcsummed ra- 
dial potential ^(u; a§, a§; j^) will play in our work the 
role played by the nonlogarithmic 5PN Padeed poten- 
tials A'^°~^°s{u; as, og; i^) (obtained by replacing a§(i^) -I- 
a!f(i^)lnu — > aQ{iy) and a§(j/) -\- a^g{iy)hiu — >■ 05(1^) in 
the formulas above) used in the previous EOB works 
[6, 22, 28]. As in those references, we shall use NR data 
to constrain, for each value of the symmetric mass ratio 
J/, the values of a§(i'), and Oq^v). To simplify this task, 
we shall take into account from the beginning the finding 
of Refs. [6, 22, 28]. The latter references found that there 
is, for each value of i/, a good EOB/NR agreement within 
a long and thin banana- like region in the (05, ag) plane. 
In view of this degeneracy between 05 and ag, we shall 
then fix the value of 05 and fit only for the (i^-dependent) 
value of ag(i^). Recent works connecting PN and/or EOB 
theory to gauge-invariant observablcs computable from 
gravitational self-force (GSF) theory have succeeded in 
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determining the j/-linear contributions to the two EOB 
potentials A{u;v) and B{u]v) [35, 38, 40, 45-47]. In 
particular, the limiting value as v Q oi the "true" 
Taylor value of a'i^{v) (defined from Eq. (7)) was found 
to be al^''^^°\Q) = 23.50190(5) in Ref. [38]. Note 
also that Ref. [47] derived some "effective" values of 
05(0) from w-global fits of the linear contribution to 
A{u,v). For example, from their best global fit, they 
obtained a§''^(0) = 23.47267. In view of these results, 
and of the fact that the parameter a§(;^) that we use in 
our EOB work is (after Fade resummation) an effective 
parametrization of the shape of the A potential, we de- 
cided to fix our a'i^{v) parameter to the simple, rounded- 
off, z^-independent, fiducial value 



c fiducial / 



v) = 23.5. 



(11) 



Finally, as the other EOB potential B{u), or equivalently 
the associated potential 



D{u; v) = A{u; v)B(u, v), 



(12) 



plays only a secondary role in the dynamics of coalescing 
binaries, and is therefore difficult to probe by using NR 
data. We used its 3FN-resummed value as obtained in 
Ref. [3], namely 



D{u; v) 



1 



1 + &vu^ + 2(23 - iv)vu^ 



(13) 



without trying to improve it by including the known loga- 
rithmic contributions appearing at 4FN and 5FN [37, 38] 
(which mix with unknown nonlogarithmic contributions). 



Summarizing: Our EOB Hamiltonian 

HEOBir, Pr,,P^p) contains only one free (i/- 
dependent) parameter, namely The 
Hamiltonian HEOBi''',Pr,-,Pip]ciQ{i')) is defined by 
Eqs. (3), (4), (10), (13), with Eqs. (7), (8), (9), and (11), 



together with pr, 
and u = 1/r. 



{A/B)^/^Pr , B{u) = D{u)/A{u) 



B. Improved EOB Waveform during inspiral and 
plunge 

Following Refs. [5, 8, 9], we describe the inspiral-plus- 
plunge multipolar waveform by the factorized structure 



insplungc 



(14) 

where we indicated the (main) arguments used in several 
factors of the waveform. Here e = 0, 1 is the parity of the 
considered multipole (i.e. the parity of ^ + m) , /i^^''^'' the 



Newtonian waveform, 5"^^ a source factor, with S^^ 

HcS or S^g = Pipl {tuiV^p) according to the parity of the 
multipole (see below for definitions). 



(0) 



(15) 



the tail factor [5, 8, 9], p^n the resummed modulus 
correction and a next-to-quasi-circular correction. 

The precise definitions of the factors entering Eq. (14), 
and of their arguments is given next. 
The Newtonian contribution reads 

= ^n'~;ic,U-)vl^'Y^--- (|,^) , (16) 

where ip is the orbital phase, ~ a, suitably defined 
azimuthal velocity, and = 
dius with ip defined as 

^(-^P.) = ^ (^) 
1 + 2j/ 



rip'^/^ a modified EOB ra- 



A 1 




(17) 



The definitions of v^p and are such that they satisfy 
Kepler's law, 1 = il^r^ = v'^ruj, during the adiabatic 

inspiral [48]. In Eq. (16), nf^^ and C£+g(i/) are numerical 
coefficients given by [5] 



(im) 



[im) 



Wni 



(18) 



'(2^+l)(^-^-2)(^2_,^j2) 



{2e + iy.\y {2i - i){e + i)i{e - 1) 



e+e-1 



-i)™xf+<^-i 



(19) 
(20) 



where X1.2 = mi^2/M . For what concerns the tail factor 

^tail 



/i^f^', Eq. (15), its main contribution, Tf™, is written as 



r(£ + 1 2ifc) ^fc 2ifc ln(2fcro) 

r(£ + i) 



(21) 



with k = ttiGHbob^, k = mfl and tq = 2GM/^. Note 
that, apart from the logarithm term ln(2A;ro), the main 
tail contribution T^,„ depends on the dimensionless ar- 
gument y = {GHeob^)^^^ : which differs from the usual 
dimensionless frequency parameter x = (GMJ7)^/'^ by 
the replacement M — ?> iJEOB- 



1. Further resummation of the residual tail phase Si,-n{y). 

The main factorized tail term Timiy) = \Tim\s''^'™ is a 
complex quantity whose modulus \Tim\ describes the tail 
amplification of the waveform modulus, and whose phase 
Tim describes the main part of the dephasing caused by 
tails. There are, however, additional dephasings caused 
by tails, which are described by the supplementary phase 
factor e^^'^"^ in Eq. (15). The residual phase correc- 
tions Si„i{y) entering the tail factor (15) were obtained in 
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Ref. [5] as a PN series in the variable y = (GHeob^)^^^- 
Here we shall use for Sgmiu) an expression that differs 
both from the one originally given in Ref. [5], and from 
its test-mass-higher-PN completion given'^ in Ref. [51]. 
More precisely: (i) we do not include the highest-order 
0(y^/^) test-mass {v = 0) PN corrections because of their 
PN-gap with respect to the last known comparable-mass 
terms; (ii) we include the 3.5PN, i^-dcpcndent, contri- 
bution to (522 (y) that can be deduced from a recent an- 
alytical computation of the PN-expanded waveform at 
3.5PN accuracy [44]; and (iii) we Pade-resum the Taylor 
series in powers of y^^^ giving Simiv)- Indeed, we found 
that the PN-expanded version of 5irn(]j) presents some 
unpleasant features (discussed below in the i ~ m = 2 
case) that are avoided if one resums 5im{y^^'^) by factor- 
izing the leading-order term and and replacing the rest 
with a suitable Padc approximant N{y^^^)/ D{y^^'^). 

Let us explain our new procedure on the (most im- 
portant) example of the £ = m = 2 phase (the others 
are listed in Appendix D). Let us start from its Taylor- 
expanded form 



7 
— f 

3' 

/ 30995 



5^^''°\y)^^y'/'-2A.y'/' + ^^ny' 



962 



V 1134 135 



7/2 



(22) 



Here we did not include the highest-order test-mass 
term (-2203/81 1712/3157r2) y^/^ that was obtained 
in Ref. [51]. On the other hand, the 3.5 PN i^-dependent 
term proportional to y''/^ is a new contribution that is 
obtained by applying the factorization of [5] to the re- 
sults of [44]. Note that this is the only genuinely new 
information given by this calculation; indeed, the real 
3.5PN contributions to /122 are already contained in the 
modulus of the EOB-resummed tail factor h^^. For the 
comparable-mass cases that are of primary concern for 
upcoming GW observations (say for v > 0.1) the 0{y'^/'^) 
contribution is numerically quite significant compared to 
the lower-order terms. To better appreciate the relative 
importance of the successive PN corrections we factorize 
Eq. (22) in a leading order (LO) part, 5\^{y) = (7/3)y3/2 

and a fractional PN-correction term, i522 = '^22''''° V'^2 ° ■ 
In terms of Vy = ^/y, the latter fractional PN-correction 
has the structure 



1 + C2W„ + C3V + C4V 



(23) 



We plot, in Fig. 1 the successive truncated PN approxi- 
mants, at IPN, 1.5PN and 2PN accuracy (i.e. up to Vy, 
Vy and Vy) for g = 1 (i^ = 1/4) and q = Q (ly — 6/49 w 
0.1224). This figure illustrates two facts: (i) the succes- 
sive PN approximants to ^22 = ^ + C2Vy + C3Vy+C4Vy + - ■ ■ 
are suspiciously different from each other; and (ii) they 



introduce rather large fractional modifications of the LO 
phase 5^^{y) = [7 /i)y^/'^ when Vy > 0.3 (which is 
reached during the late plunge). This suggests a non- 
robust behavior of the Taylor approximants in the high- 
velocity regime. In addition we have found that using 
622^^°' (y) in the generation of the FOB waveform gener- 
ates pathological features in the waveform phase in the 
very late plunge phase, compromising its accuracy of the 
phasing in a crucial region. To overcome this difficulty, 
we replace S22^^™{vy) with its (2,2) Pade approximant, 

i.e. we take P2 ['^22('yn)]- Finally, we use in defining the 
factorized FOB waveform the following resummed ver- 
sion of the 522{y) phase: 



522{y)^5\^{y)P^ 



^22 (Vy) 



7 .PQ + PlVy + P2Vy 

—v:, 

3 ^ 



Po 



PlVy+p'^vl 

(24) 



where 



The explicit expressions of the in- 



dependent Pade coefficients Po{v),pi{i^),P2{i^),p'2{^) will 
be found in Appendix D. Note that this Pade representa- 



tion degenerates as — ?> 0, and yields P| i522 {vy) — )• 1; 

this occurs because the definition of this Pade approxi- 
mant crucially depends on having a non- vanishing 3.5PN 
contribution. Figure 1 compares the Pade-resummed 
(522 (fy) to its successive Taylor approximants. This figure 
suggests that the Pade approximant represents a reason- 
able "average" of the successive Taylor approximants. 
We found that the (known) successive PN approxi- 



mants to 621^^' 



JTay: 

"33 



lor 



and (5jf , exhibited a rather 



non robust behavior similar to that of 622 ■ We there- 



Taylor 
22 

fore decided to Pade resum them, using now (1,2) Pade 
approximants, in view of the available PN knowledge. 
For the other residual phase corrections, 632, 64^, with 
m = 1, ... ,4 and S^^, there is too little PN information 
to try a resummation, so that we keep them in their un- 
resummed Taylor-expanded form. See Appendix D for 
details 



2. Further factorized corrections to the waveform: pimiv^p) 
and L^-^C^ 

Let us first emphasize that, as in our previous 
work [22], we shall use as argument in the modulus cor- 
rection pi,n (to replace the generic variable x used in 
[5]) the quantity = v^ = (r^fl)^ defined above. By 
contrast, Ref. [28] uses x = (MQ)^/^ as argument in 
the Pirn's. The pem's that enter Eq. (14) are taken at 
the complete 3"*"^ PN approximation (as done in previous 
work [26, 27, 31, 42, 52]), i.e., by completing the 3PN- 
accurate, ^/-dependent results of Ref. [5] by the = 0, 
5PN-accurate, terms obtained by Fujita and Iyer [51]^. 



^ In successive steps, this computation has been recently pushed 
to the remarkable 22 PN order [49, 50]. 



^ Let us mention in this respect that Fujita has recently extended 
the PN accuracy of the 1/ = pirn's to the 14PN level [49]. 
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where the rii's factors are chosen here to be 



0.15 

v.. 



FIG. 1. (color online) Comparing the Taylor-expanded 822 
with its (2,2) Fade approximant for two mass ratios. 



Note that in doing so we are taking into account more 
test-mass terms in the pim^s than was done in Ref. [28], 
which was stopping one PN order earlier for P33, p^i, 
/94m, and two PN-order earlier for P5m, pem and pjm- 
For completeness we list in Appendix D the explicit ex- 
pressions of the pim^s that we use. As we said, one must 
replace the generic variable x used in these expressions 



by 



Let us now discuss the structure of the final, NQC fac- 
tor h^^^ in the factorized waveform, Eq. (14), as well 
as the procedure we shall use to determine (from NR 
data) the values of the coefficients aj"^ and bj"^ entering 
this NQC correction factor We shall adopt here a 

more elaborate NQC factor h^^^' than what was consid- 
ered in previous EOB literature. In particular, for each 
multipole {i, m) this NQC factor depends on 6 real pa- 
rameters, 3 for the amplitude, of™, i = 1,...,3, and 3 
for the phase ^f™, i — 1, ... ,3 and reads 



fNQC 



cxp 



i=i 



"1 = -FT 
Vrf2/ 


(26a) 




(26b) 




(26c) 


Pr. 


(26d) 




(26e) 


"6 "SPr. ■ 


(26f) 



Here, the superscript (0) on the right-hand side of the 
definition of 712 means that the second time derivative 
of r is evaluated along the conservative dynamics (i.e. 
negecting the contributions proportional to J-, see Ap- 
pendix A for a discussion). 

For each multipole, the 6 parameters af™ and bf" en- 
tering Eq. (14) are determined from NR data by imposing 
that the EOB waveform /if°^(t^°^) (which is a func- 
tion of the EOB dynamical time i^*-*^) "osculates" the 
NR waveform hf^{t^^) (which is a function of the NR 
retarded time t^^) around some "NQC-detcrmination 
point" , which corresponds to the time p^ak on the EOB 
time axis, and to the time t^^,- on the NR time axis. In 
this work wc shall use as NQC-dctermination point on the 
EOB time axis the EOB dynamical time tf^^^^^ when the 
EOB orbital frequency ri(t^*-'^) reaches its (first) maxi- 
mum. This EOB time was often referred to, in previous 
works, as the "effective EOB light-ring crossing time", 
because, in the test-mass limit, it does correspond to 
the dynamical time when R{t^'~'^) — 3M, and, in the 
comparable-mass case, it is very close to the time when 
^j-^EOB-j ci-osses the formal EOB analog of the light ring. 
Here, to avoid confusion, we shall call it the f2-peak time, 
and denote it as t^poak- We shall discuss later the choice 
of time, say t^^^., corresponding on the NR time axis to 
the EOB time i^OB ^. 

The degree of osculation between the EOB and NR 
waveforms is defined by separately imposing a contact 
between the amplitudes, and the frequencies, of the two 
waveforms at the NQC-determination point t^p^ak ^ 
t^^j.. We do not constrain the relative phase of the EOB 
and NR waveforms. Explicitly, we impose the following 
six conditions 



4 EOB /.EOB \ 
^Im I'^Opcak; 


/iNR/.NR \ 
— ^Im l^'-extr/i 


(27a) 


4EOB/.EOB \ 
I '-a peak/ 


4NR/.NR \ 


(27b) 


4EOB/.EOB \ 


;iNR/.NR \ 
~ ^im I'-extr/J 


(27c) 


, , EOB /.EOB \ 
^im l^^Opcak; 


_ , ,NR/,NR \ 


(27d) 


, -EOB/. EOB \ 
^im l^^Opcak; 


_ •,NR/,NR \ 
— ^im I'-extrJ) 


(27e) 


,-,EOB/.EOB \ 


_ --NR/.NR N 
— ^Im I'-oxtrJ) 


(27f) 



(25) which yield two separate 3x3 linear systems to be solved 
to obtain the af™'s, and, separately, the fof^'s. 
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Note that the values of the a^'^'s affect the modulus 
of the inspiral-plus-plunge waveform, which then affects 
the computation of the radiation reaction force (through 
the angular momentum flux, see below). In turn, this 
modifies the EOB dynamics itself, and, consequently, the 
determination of the (a^™,6j™)'s. This means that one 
must bootstrap, by iteration, the determination of the 
(a^™, &^'")'s until convergence (say at the third decimal 
digit) is reached. This typically requires three iterations. 
In previous work only the dominant (2, 2) NQC cor- 
rection was included in the radiation reaction (though 
they were all taken into account when finally comparing 
EOB and NR waveforms). Here we shall follow the 
same simplifying prescription, though we have explored 
the effect of including also the subdominant (2,1) and 
(3, 3) NQC corrections to the fiux. We found that their 
effect amounts only to a small change in the NR de- 
termination of the "good values" of Og (see Appendix C). 

Summarizing: Our EOB waveform is given by 
Eq. (14) and employs the rcsummation of residual ampli- 
tudes Sgm as in Eq. (24). The NQC correction is defined 
by Eqs. (25)- (26) with constants determined from NR 
data by Eq. (27). 



C. EOB waveform during merger and ringdown 

One of the specificities of the EOB formalism is to 
construct a complete waveform, covering the full pro- 
cess from early inspiral to ringdown, passing through 
late inspiral, plunge, and merger. This is done by at- 
taching a sum of quasi-normal modes (QNM) to the end 
of the plunge waveform. The procedure for doing so 
has improved over the years [2, 8, 53]. Here, we use 
a new way of extending the inspiral-plus-plunge wave- 
form to describe the merger-plus-ringdown subsequent 
signal, which fits with the NQC-determination procedure 
we have explained above. Our new procedure for, simul- 
taneously determining NQC corrections, and attaching 
QNM's, is motivated by the findings of Bcrnuzzi, Na- 
gar and Zenginoglu [54] in the extreme mass ratio limit 
{i^ <^ 1). We shall discuss the rationale for this procedure 
in the next Section. 

The merger-plus-ringdown signal is described, for each 
multipole £m, by a sum of N QNM signals of a final Kerr 
black hole (of mass Mf and spin parameter a/), say 



n 2 \ ^-1 

I I ringdown/, N _ \ ^ r<irn ^-cr^ ' 
^ n=0 



peak) 



(28) 



where cri''^™ 



is the complex frequency 



of the nth QNM of multipolarity Hm and C„ complex 
constants. 

In this work, we use N ~ 5 positive frequency (wf™ > 
0) QNM's. These complex frequencies are functions of 
the mass Mf and spin parameter a/ of the final hole [55] . 



For Mf and a/ we adopt the fit to the numerical results 
given in Eqs. (29) of [28], 



^ = l+ (W^-l)i/ - 0.4333i^2 _ 0.4392:/^ (29) 



12 V- 3.871i^^ + 4.0281^''. (30) 

IVl 

The procedure we shall use here for matching the ring- 
down signal (28) to the inspiral-plus-plunge signal (14) 
is similar to the ones used in previous EOB work [54] 
though it differs in a significant way from the one used 
in [28]. Namely, contrary to the latter reference, the at- 
tachment (along the EOB dynamical time axis 
the QNM signal (28) to the NQC-corrected inspiral-plus- 
plunge signal (14) is done, for each multipole €m, at the 
time 



iEOB _ .KUB _ j-BUB /oi \ 

'•fmQNMattachmcnt ~ "-f mmatching ''ilpcaki \'^^> 

where we recall that t^p^ak denotes the EOB dynamical 
time where the EOB orbital frequency reaches its (first) 
maximum. Note in particular that i^p^^k does not de- 
pend on the considered multipolarity £m, so that we are 
attaching the QNM's corresponding to all the difi'erent 
multipolarities at the same EOB dynamical time. 

To complete the description of our QNM attachment 
procedure it remains to say that we determine, for each 
multipolarity im, the values of the N complex coeffi- 
cients by requiring that the (NQC-corrected) EOB 
inspiral-plus-plunge waveform /i^'^P'""*''(t^°^), Eq. (14), 
coincides with the QNM sum (28) at N points, say 
^1, ^2, • • • , tN, forming a regularly spaced "comb" on the 
^EOB g^^[g^ centered on i^pgak- Such a "matching comb" 
is specified by choosing its total length, say 



.EOB 



.EOB 



. match 



(32) 



D. 



Improved radiation reaction: including horizon 
absorption and a radial component 



Let us now turn to our improved description of the 
radiation reaction force entering the EOB dynamics. 
Note that we have included in the equations of motion (5) 
not only an azimuthal radiation reaction J-'^ (as in all 
previous EOB works), but also an explicit radial contri- 
bution Tr, ■ We have improved the analytical description 
of both components of T. Let us discuss them in turn. 

The azimuthal component, J^ip, of the radiation reac- 
tion force describes the loss of the orbital angular momen- 
tum of the system during evolution. Indeed, Hamil- 
ton's equation for p^p reads 



dt 



(33) 



where ~ T^jv. 
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Following a standard BOB practice (since Ref. [2]), we 
require that the loss of orbital angular momentum be 
balanced by the instantaneous flux of angular momen- 
tum leaving the orbital system. In previous EOB work, 
one took into account only the flux of angular momen- 
tum in the form of GW's at infinity. However, there is 
also a flux of angular momentum which is drained out of 
the two-point mass orbital system by penetrating within 
the two horizons of the moving black holes. [The lat- 
ter flux is transformed from the orbital form measured 
by to some intrinsic spin-angular momentum of the 
holes; from the point of view of the orbital p^p this repre- 
sents a loss that must be accounted for by an additional 
contribution to J^^p.] We shall include here such an addi- 
tional horizon- absorption flux by using the recent work 
of Nagar and Akcay [41]. The corresponding effect is 
rather small and, in a PN sense, starts only at the 4PN 
level [56, 57]. Reference [57], using a leading-order (New- 
tonian) approximation both to the phase evolution and 
to the horizon flux has estimated that, in the nonspinning 
case (that we consider here) , the inclusion of the horizon 
flux entails an additional dcphasing at i? w 6M smaller 
than 0.01 rad for mass ratios 1 < g < 4. On the other 
hand, recently Bernuzzi, Nagar and Zenginoglu [42], us- 
ing an EOB description of the phase evolution together 
with an improved estimate of the horizon flux (resum- 
ming higher effects), have found significantly larger de- 
phasings (accumulated over the last 20-30 orbits) than 
those estimated in [57]. Within the EOB model that we 
use here we confirmed the findings of Ref. [42] . For in- 
stance taking the most relevant case q — 6 with initial 
separation ro = 15 (corresponding to ^ 27 orbits up to 
merger, see Table II below) the effect of horizon absorp- 
tion entails a dephasing A^cp = cj)^^-^ —cj)-^ ^ 0.12 rad at 



EOB 
£1 peak 



that increases up to 0.18 rad during ringdown^. 
Such dephasings are quite significant for the EOB/NR 
comparison that we shall perform below. This is why we 
decided to include the horizon contribution to the angu- 
lar momentum flux. 

It is convenient to decompose J-^ as the product of 
the usual quadrupolar GW flux (expressed in terms of 
and of the orbital frequency fi — dip/dt) and of a 
supplementary dimensionless correction factor (of the 1 + 
0(a;)-type) : 



contribution, and can be further written as 

fix; v) = /^(x; I.) + (1 - 4;. + 2i,^)x^ f" {x; v), (35) 

where each function p'^ '^^(x\ v) is of the 1 + 0{x) type 
and is defined by dividing by the corresponding £ = m = 
2 LO contribution, namely 



22 



(36) 



Here, F^'^''^'^ is the total asymptotic {J^) and horizon 
{H) energy flux for circular orbits summed up to mul- 
tipole i = Cax, while F^'^'^ = {32/5)i^'^x^ is the LO 
(or "Newtonian") quadrupolar (asymptotic) energy flux, 
and ^2^'^° = (32/5)1/2(1 ~ 4„ + 2v^)x^ = ^^(l - + 
2j/2)F^'^ the LO quadrupolar horizon flux [56, 57]. In 
the EOB model one uses suitably factorizcd expressions 
for the multipolar fluxes _F^m ' to resum and improve 
them with respect to standard PN-expanded expressions 
in the strong-field, fast-velocity regime. In the case of 
the multipolar asymptotic fiux Fj^^^, this factorizcd fiux 
is simply defined (as first proposed in [22]) by squaring 
the corresponding factorizcd multipolar waveform of [5] , 
recalled above. An analogous procedure for the multipo- 
lar horizon fiuxes, F^ was introduced in Ref. [41] and 
compared with Regge- Wheeler- Zerilli numerically com- 
puted horizon fiuxes in Ref. [42] . [Here, we are consider- 
ing nonspinning binaries.] 

The horizon and asymptotic energy fiuxes along circu- 
lar orbits are then written as multipolar sums, say 



(x; v), (37) 



£=2 m=l 



where F, 



,H.i 
^1 



sums the two equal contribu- 



tions corresponding to +m and —to (to 7^ as the m = 
contributions vanish for circular orbits). 

Inserting in the (circular) asymptotic multipolar fiux 
contribution. 



(38) 



32 



(34) 



the factorizcd waveform (14) yields 
2 



F, 



^ ^ (39) 
Here the function f {x;^) = 1 + 0(x) (taken with the ^^^^^ p{N,e) ^^^^^^ inserting the Newtonian-order 



argument x = w^) is the reduced flux function. It can 
be defined, for a circularized binary, as the ratio between 
the total energy fiux (including the horizon flux) and the 
^ = TO = 2 asymptotic energy flux. In our case this func- 
tion is given by the sum of an asymptotic and a horizon 



waveform in (38), and where each subsequent factor is 
the squared modulus of a corresponding PN-corrcction 

factor entering (14); e.g., f''^^''^^^'^'^ 



"■em 



1 + ■ Let us mention that \Tim{y)\'^ can 

be explicitly written in the simple form 



^ Note that one has <f> ~ 1.6 X 10 ''at the initial separation 
ro = 15, which is neghgible compared to the dephasing accumu- 



lated during the subsequent evolution. 



\Tem{y)\ 



Airk 



ii\f 1 



^ s—1 



n i-^' + (2^)' 



(40) 
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TABLE I. Coefficients of our hybrid 1+^PN-accurate 
pf^ix; v) functions as given by Eq. (43). 



I 


m 


Im 
Cl 


C2 


tm. 


4™ 


2 


2 


4-211^+271/^ -Si'^ 
4(l-4i/+2i/^) 


4.78752 


26.760136 


43.861478 


2 


1 


0.58121 


1.01059 


7.955729 


1.650228 



Similarly the horizon partial multipolar fluxes are writ- 
ten in factorized form [41] 

(41) 

where pfj^ix; z^) = 1 + 0{x) are the residual amplitude 
corrections to the horizon waveform. Following Refs. [41, 
42] we use a 1+^ PN approximation for pf^{x\ v) and we 
include only the £ = 2 contribution in Eq. (41) (i.e., we 
fix Cax = 2 in Eq. (37)). 

Finally, this means that the fractional horizon cor- 
rection (before multiplication by the additional factor 
1 - 4i^ + 2i/2) in Eq. (35) is of the form 



x^f"{x- v) 



+x° 



where S^^ = Hcs, S^s = P^l {'^uiV^), and where we use 
4PN accurate expressions for p|f^(a;; v), 

pf^ix; iy) = l + c^^x + S^x" + c'^x^ + (43) 

with values for the needed 1 — 2 coefficients cf™, i = 
1, . . . ,4 hsted in Table I. 

Let us finally come to discussing the radial compo- 
nent Tr, of the radiation reaction force. Such a contri- 
bution was generally neglected in previous FOB papers, 
or replaced (e.g. in [10, 28]) by an expression which was 
not consistently derived. Recently, Bini and Damour [43] 
(building on previous work by Iyer and collaborators [58- 
60] ) have shown that consistency with the usual FOB def- 
inition of T^p (as being equal to minus the instantaneous 
flux of angular momentum) required a specific form for 
J-r, which differed from previously used expressions. 

The result of Rcf. [43] that we use here has the form 

Tr. = -^—J'A^ + ciiiy)u + c2{yV) (44) 

where the coefficients entering the 2PN correction 
read [43] 

227 1957 

ciiv) = v^ , (45) 

^ ' 140 1680' ^ ' 

, , 753 2 165703 25672541 

C2{v) = 1' H v . (46) 

^ ' 560 70560 5080320 ^ ' 



E. Post-post-circular initial data 

The construction of initial data for the FOB dynam- 
ics has been refined in a series of works [2, 9, 10, 20]. 
Here we shall use the post-post-circular prescription, in- 
troduced in 2007 (see Sec. Ill B of [9]), and then used in 
all subsequent FOB-related works by our group [15, 20- 
22, 24-26, 31, 42, 52]. This choice allows one to start 
the FOB dynamics (with negligible initial eccentricity) 
at a frequency that is compatible with the initial fre- 
quency of the NR waveforms we shall use here (MW22 ~ 
0.0345 approximately corresponding to initial separation 
i?o ~ 15M, see Table II below). Note that, by contrast. 
Pan et al. [28], who use the less accurate post-circular 
initial data of RcL [10], start their FOB runs at an ini- 
tial radius Rq > 50M (corresponding to an initial GW 
frequency M0J22 < 0.005) in order to get a good circular- 
ization of the dynamics at the frequency where numerical 
simulations start. 

For completeness, let us review here the construction 
of post-post-circular initial data for a given relative ini- 
tial separation tq. We introduce a formal book-keeping 
parameter e (to be set to 1 at the end) in front of the 
radiation reaction T^p in the FOB equations of motion. 
The quasi-circular inspiralling solution of the FOB equa- 
2tions of motion can then be formally expanded in powers 
of e as 



Pr, 



'1 



^fc2(r)+0(£^)) 



r^i(r) + 0(£3). 



(47) 
(48) 



Here, jo('') is the usual circular approximation to the in- 
spiralling squared angular momentum as explicitly given 

by 



Jo(0 = - 



A'{u) 
u^A{u)'^ 



(49) 



where the prime means d/du (recall u = 1/r). The or- 
der e-approximation to p^, j i-e. 7ri(r) ( "post-circular" ) is 
then obtained by approximating the left-hand side (l.h.s) 
of Fq. (6c) by dp^/dt w djo{r)/dt = {djo{r)/dr){dr/dt). 
This determines dr/dt and thereby a corresponding value 
of Pr, using Eq. (6b) (where we neglect the p^, contribu- 
tion). This leads to the following explicit expression for 
71-1 (r): 



STTi{r) 



1/2 



dr 



(50) 



J 



where the subscript indicates that the r.li.s. is evalu- 
ated at the leading circular approximation e — >■ 0. The 
post-post-czrcw/ar approximation to (term e^fc2 above) 
is then obtained by approximating the l.h.s. of Fq. (6d) 
by 



dpr, d'Ki{r) dr 
dt dr dt ' 



(51) 
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TABLE II. Post-post-circular initial data for EOB dynamics 
that we shall consider in this paper to assure negligible initial 
eccentricity. They are obtained with the choice a% = 23.5 and 
al{u) = -122 147(1 - Av). 



q 


V 


ro 




Vr 




1 


0.25 


16 


4.42468388 


-0.00101207 


-0.00088971 


2 


0.2 


15 


4.31684806 


-0.00113064 


-0.00098465 


3 


0.1875 


15 


4.31889797 


-0.00096446 


-0.00083930 


4 


0.1600 


15 


4.32052780 


-0.00083019 


-0.00072202 


6 


0.1224 


15 


4.32276240 


-0.00064296 


-0.00055874 



where the radial derivative d'Kx (r) / dr is numerically com- 
puted. This transforms Eq. (6d) in a linear equation 
for p^, which leads to an explicit expression for the 
dependent correction e^fc2(r) introduced above. In solv- 
ing for pfp we keep for additional accuracy the contribu- 
tion proportional to p^^ ~ e^T:\ (r) . 

Table II lists the post-post-circular data (as a function 
of rg) obtained by this procedure, as we have used them 
in the present study. Note that these values mainly de- 
pend on the parameters entering the A function, (a§, a§), 
and depend almost negligibly on the values of the NQC 
parameters a|™ entering J-^p that appears on the r.h.s. 
of Eq. (50). Actually, the values listed in Table II were 
computed by keeping only the {a^,a^) NQC contribu- 
tions. 



F. Analytically unknown parameters, choices to be 
made, NR-completion of the EOB model 

Let us summarize the parameters entering the con- 
struction of our EOB model, emphasizing which param- 
eters contain important dynamical information, which 
ones are already known with sufficient accuracy, which 
ones depend on reasonable choices we can make, and how 
NR data can be used to complete the EOB model by de- 
termining the various parameters. 

At face value, the EOB model defined above depends 
on quite a few analytically unknown parameters, namely: 
a'i(y), a§(i^), the 6 NQC parameters (af™,&f™) for each 
waveform multipole, the values of the mass Mf and spin 
a/ of the final black hole, the number N of QNM modes 
used in the ringdown signal, and the width A'"''*'^'^ of the 
QNM matching comb. 

Our attitude towards the use of NR data to complete 
the EOB model by determining these parameters is the 
following: 

(i) As already said, we think (in view of previous EOB 
results [22, 28]) that it is a reasonable choice to 
impose some a priori relation between a§(i/) and 
a%{v), so as to look only for one free dynamical pa- 
rameter. Here we shall fix a^iv) to the simple value 



a|(i/) = 23.5, Eq. (11). This leaves only a§(i/) as 
free parameter. We shall discuss below (see Sec- 
tion V) how the nonperturbative information con- 
tained in NR phasing data can be used to determine 
the value of a%{v)^ in a way which is nearly decorre- 
lated from the uncertainties in the determination of 
other parameters. Let us already indicate here the 
very simple final result we shall get for this EOB 
parameter 

ag(i^) = -122 + 147(1 -4i/) (52) 
from Caltech-Corncll-CITA data 

We think that the NR determination of ag(i^) leads 
to important information about the conservative 
dynamics of binary black holes (as we shall illus- 
trate below); 

(n) Concerning the NQC parameters (af"*,&f"*), the 
procedure explained above reduces their determi- 
nation from nonperturbative NR data to a single 
choice, namely that of the time t^^-^ on the NR (re- 
tarded) time axis corresponding to the EOB time 
^Opcak (which can be thought of as defining the 
"EOB merger time"). The choice of on the NR 
time axis is not a matter of convention, but has (a 
priori) important physical consequences. It must 
be done by combining information coming both 
from comparable-mass NR simulations, and from 
extreme-mass-ratio ones. For reasons that shall be 
discussed below, we shall choose, for each mass ra- 
tio V, a specific value of t^^^{v) given by 

t^lil^) = t^A^. peak(^) + /(^) (C pcak(^) " ^A^, pcak(^)) 

(53) 

where 

./M = ^(l + 7(l-4i.)), (54) 

and where t^^pcaki^) ^he NR time when the 
NR quadrupolar amplitude reaches its peak, and 
t^^^-pcakiy) the NR time when the quadrupolar fre- 
quency has an infiection point. Here f{i') varies 
between /(O) = | and /(i) = as v varies be- 
tween and i. t^^^{v) always lie on the right 

of (i.e. later than) the NR time i^f^ poak('^)- We 
shall extract nonperturbative information from NR 
data by computing from the various multipolar NR 
waveforms a certain number of derivatives of their 
amplitudes and frequencies at the extraction point 

(iii) Building on previous work, we shall use the simple 
(NR-based) analytical fits (29) for the mass and 
spin of the final black hole. Note, however, that, in 
principle, the EOB model (when NR-completed by 
NQC corrections up to merger) does yield, by itself, 
predictions for Mf and af [2, 61]. This might be 
useful in cases (e.g. with large, preccssing spins) 
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where one does not have in hands accurate ana- 
lytical fits for the characteristics of the final black 
hole. 

(iv) We shall use here = 5 QNM's, and as ex- 
plained below, we shall fix A™^"='^ = 0.7M for all 
multipoles. Note, that by contrast, Ref. [28] uses 
iV = 8 QNM's and introduces "pseudo QNM's", 
and much larger matching intervals arc employed, 
which also vary with £m. [E.g., the latter reference 
uses 5A/ and A^^ = 12M.] 

(v) Let us finally note that (contrary to [28]) we shall 
not introduce adjustable parameters in the wave- 
forms, nor shall we introduce special modifications 
to improve the behavior of some subdominant mul- 
tipoles. 

III. NUMERICAL RELATIVITY 
INFORMATION, AND DIAGNOSTIC 

A. Overview of numerical waveforms data 

The NR data we use here to complete the EOB wave- 
form were obtained with the Spectral Einstein Code 
(SpEC) developed by the Caltech-Cornell-CITA collab- 
oration [62-66]. Specifically, we used the waveforms 
recently published in Ref. [34], coming from simula- 
tions of nonspinning black hole binaries with mass ra- 
tios q = m2/mi = (1,2,3,4,6). Before their publica- 
tion, these data were already used in some EOB /NR and 
PN/NR comparisons [28, 33, 67]. We address the reader 
to Ref. [34] for all technical details about the numerical 
setup and estimates of the accuracy. Here we only re- 
call that these are the longest published waveforms to 
date (together with the 33 orbits, equal-mass waveform 
of Ref. [67]), with a number of gravitational wave cycles 
up to merger (here conventionally defined as the maxi- 
mum of the modulus of the quadrupolar metric waveform 
\h^f-\) respectively A^gw = {33,31,31,31,43}. We made 
use of two different types of waveform data: curvature, 
-01^, and metric, hg„i, extrapolated to infinite extrac- 
tion radius. Indeed, the metric waveform him was also 
directly extracted from the numerical spacetime using a 
Regge-Wheeler-Zerilli- based (RWZ) approach ^, sec Ap- 
pendix of Ref. [23] for a discussion. 

B. Estimating the NR Qu}{uj) function for the 
curvature waveform 

In this subsection we shall explain how we extracted 
from NR data a useful, intrinsic measure of the NR phase 



This type of RWZ approach was initiated by Abrahms and 
Price [68] and first implemented in the form of Ref. [23] in 
Refs. [69, 70]. 




Mu [curvature] 

FIG. 2. (color online) Top panel: raw NR, curvature wave- 
form, data; smoothed data and fit. Bottom panel: difference 
between smoothed data and the fit. 
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FIG. 3. (color online) Fitting the rescaled Qu,, curvature 
waveform, function. The top panel contrasts the smoothed 
data with the outcome of the fit. The bottom panel shows 
their difference. 



evolution, namely the Qu>{oj) function. This function is 
a convenient version of the "intrinsic phase acceleration" 
function a{Lo) introduced in Ref. [9], which was defined 
such that dbj/dt = a{u)). This function is an "intrinsic" 
measure of the time-domain phase evolution in the sense 
that it is independent of the two shift ambiguities that 



12 



TABLE III. Coefficients entering the fitting function for tlie rescaled, curvature waveform, Q^}, Eq. (62). In the second column 
we also report the frequency interval A/(aji,aj2) on which the fit was performed. 



g 




MUJ2 


m 


n2 


na 


n4 


ns 


di 


d2 


da 


1 


0.03877 


0.29654 


-27.88757 


256.94609 


-1053.85269 


1926.40123 


-1274.57280 


-6.60927 


47.87468 


-104.35366 


2 


0.04133 


0.29709 


15.51565 


-372.20973 


1725.17714 


-3145.40474 


2105.30901 


-15.77371 


90.80420 


-103.95952 


3 


0.04476 


0.29642 


6.50413 


-243.11043 


1108.15054 


-1913.96522 


1193.58571 


-14.56312 


79.22950 


-111.26577 


4 


0.04819 


0.29671 


0.52391 


-172.68858 


806.19352 


-1350.57200 


797.20936 


-14.39733 


78.72314 


-124.83300 


6 


0.04280 


0.29720 


7.18353 


-247.53679 


1096.28420 


-1833.39721 


1090.77236 


-14.59256 


75.97063 


-113.64331 



affect any time-domain phase, 4'{t): an arbitrary phase 
shift (/) — > 4> + c, and an arbitrary time shift t ^ t + t. 
The Quji^^) function is defined as 



Note that this definition is equivalent to saying that the 
time-domain phase accumulated in the frequency interval 
(wi,W2) is given by the integral 

't>{uji,uj2) = / Qu^dlwuj. (56) 

The function Quj{oj) has proven to be a very useful 
diagnostic of phase evolution in recent EOB/NR com- 
parisons of binary neutron stars [24, 25, 27]. Note that, 
in the definition, uj can be the frequency either of the 
curvature waveform, or of the metric one (thereby defin- 
ing two different, though numerically close, functions). 
In general, one only considers the frequency of the dom- 
inant quadrupolar waveform, though one can also study 
the Qu:{(^) function of any {i,m) multipolc. Note also 
that we are here considering the phase acceleration of 
a time-domain phase. One can also usefully consider 
the frequency-domain counterpart of (5lj(w), defined as 
(5cj(a;) = uj^(Pip{uj)/doj'^, where tp{uj) denotes the phase 
of the Fourier-transformed waveform. In the station- 
ary phase approximation, Q^^{uj) is simply equal to the 
time-domain Qcj(a;) (see, e.g., Eq. (17) in [31]). 

Let us now discuss how to accurately estimate Qu:{'-^) 
from the numerical data, in spite of the loss of accuracy 
associated to the fact that its definition (55) involves the 
computation of two derivatives of the phase 4){t). We 
consider the -012 curvature waveform extrapolated to in- 
finite extraction radius, decompose it in amplitude and 
phase with the convention 

^22 = 1^-1216-'*^= , (57) 

and consider as frequency in the definition (55) the cur- 
vature quadrupolar frequency: w = 022- 

It is somewhat of a challenge to get an accurate 
out from numerical data. For example, in the case of 
binary neutron star waveforms, Refs. [24, 25, 27] argued 
that the successive straightforward differentiation (using 
finite-differencing, 4th-ordcr stencils) of the numerical 



data is unable to get this information correctly, so that a 
suitable fitting of the GW phase was necessary to obtain 
something qualitatively and quantitatively correct. For 
general binary black hole simulations, due to the much 
higher resolution involved as well as due to the highest 
finite differencing operators used, direct differentiation 
could be more meaningful than in the binary neutron 
star case. This should be even more true for SpEC data, 
since they are expected to be particularly accurate. 

Therefore, as a first step we directly computed Q^j 
from the raw data simply by finite-differencing (j) twice 
to get UJ and to, i.e., applyng twice a Ist-derivative finite- 
differencing operator with 4th-order stencil. The result 
of this first step is shown, for q = 1 data, as a dashed, 
light-gray line in Fig. 2 (see also the close up). The figure 
shows the presence of high frequency noise which pre- 
vents one from using this diagnostics as is for reliable 
quantitative estimates. 

To improve on this, and get a quantitatively useful 
estimate of the Qu) curve, we applied three more steps. 
First, in order to eliminate the high-frequency noise, we 
smoothed w(t) with a Sgolay filter. Second, we com- 
puted the time derivative of the smoothed w(i), and 
then smoothed again that derivative with a Sgolay fil- 
ter. These two steps succeeded in strongly reducing the 
high-frequency noise in the curve (thick line in Fig. 2, 
blue online). However, there remained a low- frequency 
residual oscillation in the resulting Qi^curve (evident in 
the inset of Fig. 2). We do not know the precise origin 
of this residual oscillation (it might either be related to 
some small residual eccentricity in the waveform or con- 
nected to the extrapolation procedure), but we think it 
is of spurious numerical origin and that it docs not have 
any actual physical content (note that such an oscillation 
is not present in the EOB Q^^ curve). 

This led us to our third step: a fitting procedure of the 
Quj{^) function. To implement such a fitting procedure, 
it is convenient to first normalize the curve with re- 
spect to its leading-order, Newtonian part, 

Q^(c.) = A2-V3^-5/3^ (58) 

thereby factoring out the blowing up of Quj{'^) at low 
frequency. The normalized function 

Qu,{oj)^Q^/Q^ (59) 
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stays of order unity on the full frequency range (and 
— >■ 1 for w — >■ 0) and is a better starting point for 
any fitting procedure. Then we use as fitting template 
for a general analytical structure consistent with the 
structure of predicted by PN theory in the adiabatic 
approximation. More precisely, the 3.5PN-accurate ex- 
pansion of is a Taylor expansion in half-integer powers 
of X = (Mrj)2/3 (modulo some logarithmic corrections) 
that reads 

(60) 

where 

&2 
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This motivated us to fit the smoothed version (coming 
out of the first two steps) of the numerically computed 
Qi^{ll!) with a Fade- type function of the form 



1 + niXuj + n2xl!'^ + n^xf^ + n^x'^J'' + n^x_ 



5/2 



1 + diXu. 



(62) 

where x„ = (A/w/2)2/3. 

Let us now illustrate the result of performing this 
three-step evaluation of the numerical Qu]{uj) function. 
The top panel of Fig. 3 shows, for g = 1, the three 
successive estimates of the numerical Q^j'- the raw one 
(dashed line, featuring many large spikes) , the smoothed 
one (solid line), and finally the fit obtained using the 
template (62). Note that all those curves are plotted 
versus Muj. The bottom panel of the same figure shows 
the difference AQ^(w) = Q^™°°*'^<"i(a;) - Q5'(a;) between 
the smoothed data and the fit. Note that this differ- 
ence is oscillating around 0, which indicates that the fit 
has been effective in averaging away the low-frequency 
oscillation remaining after having smoothed the high- 
frequency noise. The procedure works in the same way 
for the other mass ratios, and for each one the difference 
AQ(^(cj) nicely oscillates around zero. 

We list in in Table III, for all mass ratios, the fitting co- 
efficients of the smoothed numerical to the template 
Eq. (62). Note that this list of coefficients provides a 




FIG. 4. (color online) Hierarchy of important points of 
the test-mass (Zerilli-normalized) quadrupolar metric wave- 
form (divided by u), '^22/1' = {R/M)h22/ (i^V^i) around the 
merger point. The orbital frequency Q peaks at approxi- 
mately 2/3 of the time interval between the peak of the metric 
amplitude and the inflection point of the GW frequency, i.e. 
the first peak of 1022- 



TABLE IV. Time intervals t^f^ p^^i, 
ical data considered in this paper. 



NR 

^22 peak 



for aU numer- 



,NR 

'■cl>22 peak 



NR 

.422 peak 



1 

2 
3 
4 
6 

00 



0.25 

0.2 
0.1875 
0.1600 
0.1224 




3.2493 
3.4426 
3.3261 
3.5714 
3.5681 
3.8158 



convenient way of condensing the information contained 
in the NR phasing during most of the inspiral and plunge 
(indeed, our fit worked well up to frequency Mtu ~ 0.3, 
which is quite close to the merger). This packaging of 
the NR phasing information might be useful for many 
purposes, e.g., comparing various numerical simulation, 
computing the Fourier transform in the stationary-phase 
approximation, etc. 
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FIG. 5. (color online) Test-mass waveform: comparison between RWZ waveform extracted at J"^ and EOB waveform completed 
by the 6-parameter NQC correction factor to the waveform, Eq. (25). Top panels: £ = m = 2 multipole, modulus and frequency 
(left) and phasing (right). Bottom panels: t = 2 and m= \ and I = m = ?, frequency and modulus. The ringdown is modeled 
using 5 (positive-frequency only) QNMs. 



IV. REVISITING TEST-MASS LIMIT RESULTS 

A. The new information acquired from 
test-particle computations 

Before dealing with the Caltech-Cornell-CITA 
comparable-mass waveforms, we shall revisit in this 
Section the test-mass limit case v ^ \ both to motivate 
our introduction of an NR extraction point t^'Jj. differing 
from the peak of the waveforms, and to test the perfor- 
mance of the basis of functions n^'s that we shall use in 
our NQC correction factor, (25). 

State-of-the-art computations of multipolar Regge- 
Wheeler-Zerilli (RWZ) waveforms for the plunge and 
merger of a test particle (of mass /i), moving in a 
Schwarzschild background (of mass M), and submit- 
ted to a leading-order EOB resummed radiation-reaction 



force, have been presented in a recent series of works [42, 
54, 71, 72]. These works have used a recently devel- 
oped method [73-75] allowing one to combine an accu- 
rate treatment of the particle motion in the strong field 
region, with the extraction of the waveforms directly at 
null infinity {J). The findings of Ref. [54] that will be 
of direct interest for our present study are: 

(i) The extraction of the waveforms at J' allows one 
to relate the retarded time used as argument of 
the waveforms to the EOB time used in the 

dynamics of the particle (namely, one has simply 
^NR ^EOB-j^ This allows one to connect without 
ambiguity features in the waveform (such as, say, 
a peak in the modulus of /i22(i^^)) with features 
in the dynamics (such as, say, the location along 
the t^^^ axis of the maximum of the orbital fre- 
quency r2(t^'-'^)). Such a possibility is not available 
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FIG. 6. (color online)Test-mass limit: comparison between 
£ — 4, m = 1 EOB and RWZ modulus and frequency 



in comparable-mass NR simulations, because they 
do not track the light cones emitted by the center 
of mass of the binary system. In addition, even 
if they did, this would not allow one to relate the 
dynamical EOB time t^*-*^ to the waveform time 
t^^, because we would not know the exact relation 
between t^"-"^ and the NR coordinate time relevant 
for the NR dynamics. 

(ii) Using the connection between the waveform time 

and the dynamical time t^OB offered by (i), it 
was found that the waveform amplitude A22 peaks 
approximately ~ 2.56M earlier than the orbital fre- 
quency n, i.e. tEOB ^ ^ ^nr ^^^^^+2MM . This is a 
new information which conflicts with the standard 
simplifying EOB assumption of a coincidence be- 
tween the peaks of A22 and of il. The existence of 
a difference between t^^^^i, and i^^s peak was later 
confirmed in Ref. [76] and extended to the case of 
a spinning central black hole. 

(iii) Using this new information, Ref. [54] suggested to 
incorporate it in a new prescription for the deter- 
mination of the EOB NQC correction factor based 
on extracting numerical data at the NR point i^^,. 
corresponding to tf^°p^^^^, rather than^ at ^^^f.peak- 
They implemented such a prescription by imposing 
a contact at t^^^ o <Hpeak both (for the first 
time) between the modulus and the frequency of 
the waveform. They then showed that such a pro- 
cedure produced NQC-corrected EOB waveforms 



Note that, by contrast, Ref. [29] has chosen to keep, for £ = m = 
2 the NR extraction point at t^^^ peak ^'^'^ *° map it to an EOB 
time earlier than t^pg^i;. 



wich had an excellent agreement with the numeri- 
cal RWZ waveforms up to merger. 

The procedure we indicated in Eqs. (27a)- (27f) above 
is a generalization of this prescription to a contact 
requirement. We shall test below the increased accu- 
racy brought by using such a contact requirement, 
involving six NQC-parameters, instead of the contact 
requirement used in Ref. [54], which involved only four 
NQC parameters. This test will also probe the new basis 
of NQC correction functions n^'s used in Eq. (25). 



B. Zooming on the structure of the test-mass 
waveform near merger 

Before doing the latter test, let us display the finding 
(ii) of Ref. [54] by investigating in detail the structure of 
the £ = m ~ 2 RWZ waveform around the peak of the 
modulus, with the idea that a similar structure might 
hold in the comparable mass case. 

Figure 4 shows together (as functions of the wave- 
form retarded time u, which can be identified with 
the EOB dynamical time): the waveform modulus 
^22/^^; the orbital frequency fl; and the derivative 
of the GW frequency W22- Here, A22 is the modu- 
lus of the Zerilli-normalized quadrupolar metric test- 
mass waveform, ^22 = (R/M)h22/^/2A. [For a gen- 
eral multipole the Zerilli normalized metric waveform is 
= {R/M)ht^/^{t + 2){£ + l){l){£ - 1).] The fig- 
ure clearly illustrates how the orbital frequency peaks at 
a time i^plfak that is between the locations of the max- 
ima of A22 and a;22, i.e. we have the relation i^^pgak 
^H?oak < ^J?2^2 pcak- Quantitatively, given that we have 
'npcak ''A22 peak — ^•ooooooJw auQ t^^^ p^^,^ r^^^ p^^]^ — 
3.815784M we have that 



t 



EOB 



NR 



Q peak A22 peak 
TNR _ iNR 

6j22 peak A22 peak 



2.565388 
3.815784 



0.6723096 ; 



(63) 



The comparable-mass NR simulations show that the or- 
dering pj,g^j, < p^j^j^ remains true for all values of v 
(for nonspinning binaries). By continuity, one then also 
expects that the EOB orbital frequency will continue to 
peak between these two points for any value of v. In other 
words, one expects that the correspondence between the 
EOB and NR time axes should be such that the EOB 
dynamical time inpeak(^) corresponds to an NR wave- 
form time f^R,(z.) such that t^'^^^^.^^M < t^^,,{>^) < 

^dj^ pcak('^) ^'-'^ ^- convenient to rewrite these 

inequalities as 



*oxtr('^) ~ ^!A22peak(^) 



NR 



LU22 peak 



NR 



''A22peak('^)) 

(64) 

where /(v) is an unknown function satisfying the condi- 
tion that /(O) = 2/3, and expected to remain positive for 
any i/. We shall discuss our choice for the function f{i/) 
in the following Section. 
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C. Testing the improvements brought by requiring 
a contact when using the NQC factor Eq. (25) 



Ref. [54] was able to build a rather satisfactory EOB 
waveform modulus and frequency up to merger for the 
£ = m = 2 mode (and in general for all i = m modes) 
by using four NQC parameters (two for the amplitude 
and two for the phase). However, their results for the 
modulus were much less satisfactory for the other (£ ^ 
to) subdominant multipoles, such as the ^ = 2, to = 1 
one. Let us show here how the use of the new NQC 
factor, Eq. (25) (which contains six NQC parameters, 
and uses different choices for the NQC functions and 
714 ) improves the closeness of the EOB waveform to the 
numerical (RWZ) one. To be consistent with Ref. [54], 
the EOB dynamics used for this comparison is slightly 
different to the one we discussed above. Namely: (i) 
we set to zero J'^ , i.e. the horizon-absorption part of 
the radiation reaction; (ii) we also set = 0; (iii) in 
addition, the residual phase corrections Sim for 1^ = are 
considered in their Taylor-expanded form and all terms 
(up to 4.5PN accuracy) are included (see Appendix D). 



The improved EOB waveform obtained by using the 
new six-parameter NQC factor is illustrated in Fig. 5. 
The top panels refer to the ^ = m = 2 mode: frequency 
and modulus (left) and phasing (right). The bottom pan- 
els compare EOB and RWZ frequency and modulus for 
£ = 2, TO = 1 (left) and £ = TO = 3 (right). For all 
waveforms the QNM matching comb has a total width 
A = 0.7M and we use five, positive frequency, QNMs. 
The restriction to positive frequency QNM's is the reason 
why one cannot reproduce the oscillations during ring- 
down in the £ = 2, m = I mode. The improvement with 
respect to Fig. 3 of Ref. [54] is evident. Notably, the £ = 2 
TO = 1 modulus comes out extremely well (modulo the 
absence of negative-frequency modes to model the ring- 
down). The £ = m = 2 phasing remains good also dur- 
ing merger and ringdown -0.05 < A0^°^"^^ < -1-0.05 
(while the QNM matching of Ref. [54] led to significantly 
larger dephasings during ringdown). Note on the top 
right panel of Fig. 5 the behavior of the phase difference: 
it dips just before merger down to —0.04 rad, and then 
jumps up to -1-0.06 rad during ringdown. Such a behav- 
ior is a useful compromise to keep, on average, a good 
phasing through inspiral, plunge, merger and ringdown. 



Finally, to prove the robustness of the NQC determi- 
nation procedure and the accuracy of the EOB waveform 
for higher multipoles, we show in Fig. 6 the £ = 4, m = 1 
frequency and modulus. The agreement between EOB 
and RWZ waveform is again very good, modulo the ab- 
sence of negative modes in the ringdown modclization. 



V. COMPARABLE MASS CASE: ag(i/), Sr(i^), 
AND PHASING PERFORMANCE 



A. Iterative procedure for determining f^^r(i^) and 
af,{i'): overview 



After having tested the performance of the NQC fac- 
tor (25) in the test-mass limit, we now move to the 
comparable-mass case. Let us explain how we distilled 
crucial nonperturbative information out of the Caltech- 
Cornell-CITA waveform data. Our aim was to determine 
good values of the 5PN parameter 05(1/), and of the NR 
time t^^,{i^) corresponding to the EOB time ^opcak- We 
recall that t^^j.{i') is parametrized by a function /(t^), 
according to Eq. (64). Actually, the determinations of 
ag(i^), and of t^^j.{iy) are correlated, and must be done es- 
sentially simultaneously. From a practical point of view, 
we used an iterative, trial and error method. 

First, for a given mass ratio ly, and a given choice of 
NR extraction time t^^^. (chosen around merger), we ex- 
tract, from the behavior of the waveform in the immedi- 
ate vicinity of the retarded time, t^^,., a collection of NR 
waveform quantities {A^^,Af^,Af^,ujf^,ujf^,Qf^). 
[As mentioned above, these quantities are then used, for 
any given value of aQ^iy), to determine the parameters 
(af™, 6^™) entering the EOB NQC factor; i.e. the last fac- 
tor in the pre-merger EOB waveform (14).] Second, we 
study how the phase difference A^^'^^'"'^ between the so 
determined NQC-corrected EOB waveform and the NR 
waveform evolves (either as a function of frequency, or of 
time) from the beginning of the simulation up to t^^^- 
The evolution of the phase difference A0^*-'^^^ depends 
(after having chosen t^^^., and having implemented the 
previous step) only on the 5PN (i/-dependent) parame- 
ter ag(i^). We then search (for each whether there 
exist values of a'g{iy) which entail that A(/)^°^^" (ag(i^)) 
remains within the numerical uncertainty during the full 
simulation (up to t^^^)- If such a tuning of a§(i^) does 
not seem to lead to a a satisfactorily small phase discrep- 
ancy during the whole evolution, we try another value of 
the NR extraction time and repeat the two steps above, 
until we end up with a better pair (t^^, a§(i/)). 

When completed (by iteration), the above two steps 
completely define an NR-completed EOB model up to 
merger. The EOB waveform is then extended through 
merger and ringdown by attaching QNM's at the end 
of the inspiral-plus-plunge waveform, i.e. at the EOB 
time ^Qpeak (which corresponds to the NR time t^^^.)- 
This extension does not require the extraction of further 
NR information, but only requires to choose, by trial 
and error, reasonably good values of the number of QNM 
modes N, and of the total width of the matching comb 
^match around t^poak- already said, we use N = 5 
and A™^"=^> = 0.7M. 
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TABLE V. Next-to-quasi-circular (af", bf'") coefficients needed to complete the EOB for the five mass ratios considered. They 
are obtained by imposing conditions to the waveform amplitude and frequency around the merger. 



g 


£ m 




"2 


Im, 






6|™ 


1 


2 2 

3 2 


-0.0895 
0.0461 


1.7628 
2.4237 


0.01234 
-0.3320 


0.06741 
0.03990 


-1.0814 
0.9154 


-3.226 
-3.716 


2 


2 1 

2 2 

3 2 
3 3 


-0.0687 
-0.0659 
-0.0789 
-0.0146 


0.5265 
1.7504 
2.7396 
2.2059 


0.3144 
-0.0333 
-0.1671 
-0.1513 


0.2582 
0.0902 
0.0707 
0.2306 


0.7434 
-0.9977 
1.0058 
-1.3782 


0.969 
-2.223 
-2.410 
-2.995 


3 


2 1 

2 2 

3 2 
3 3 


-0.0569 
-0.0417 
-0.1267 
0.0080 


0.3208 
1.6715 
2.6530 
2.0324 


0.3761 
-0.0486 
-0.0782 
-0.0882 


0.2622 
0.1113 
0.1606 
0.2506 


0.8310 
-0.8073 
1.3992 
-1.2573 


1.6907 
-1.3840 
0.1241 
-1.0746 


4 


2 1 

2 2 

3 2 
3 3 


-0.0459 
-0.0320 
-0.1264 
0.0152 


0.1559 
1.5705 
2.3349 
1.8421 


0.4289 
-0.0186 
0.0247 
0.0095 


0.2766 
0.1284 
0.2697 
0.2724 


1.0053 
-0.6126 
1.8795 
-1.0269 


1.6771 
-1.1163 
0.0117 
-0.4681 


6 


2 1 

2 2 

3 2 
3 3 


-0.0317 
-0.0180 
-0.1054 
0.0340 


-0.0360 
1.4267 
1.7771 
1.5985 


0.5109 
0.0255 
0.2232 
0.1214 


0.2674 
0.1439 
0.3184 
0.2922 


1.1939 
-0.5001 
2.0994 
-0.8706 


2.2065 
-0.8488 
1.3339 
-0.1093 



B. Determining t^^r{v) 

We started by applying this iterative procedure to the 
equal-mass case q = 1 (i.e. v = 0.25). After trial and 
error, we concluded that, for q = 1, the coefHcient /(i^) 
in Eq. (64) could be taken to have a small value, such as 
/(0.25) = 1/12. In other words, when g = 1, t^^^ can be 
taken to be very close to the peak of the A22 modulus, as 
was indeed assumed in all previous EOB works. [Though 
we found that taking /(0.25) = works roughly as well 
as taking /(0.25) = 1/12, we still saw an advantage in 
having a non-zero, positive value of /(0.25).] By contrast, 
when considering larger mass ratios, we found more and 
more advantageous to increase the value of f{v), up to 
values of order of the test-mass value discussed above, 
/(O) = 2/3, for large mass ratios. Then, as a simplifying 
choice, we decided to assume for the v dependence of 
/(j^) a simple linear behavior between the two extreme 
values for 1/ = and v = 0.25, in the form 

/(j.) = /(0.25) + - /(0.25)) (1 - 41.); (65) 

which yields, when using /(0.25) = 1/12 and /(O) = 2/3, 
the explicit expression 

/(^) = ^^- (66) 

Having so chosen t^^^liy), we measure, for each 
{£,m), on the NR mutipolar waveform the vector 



at e,.(-)- Then, for 
any value of Og, we first compute the EOB dynamics, 
then we solve the linear system given by Eqs. (27a)-(27f) 
to obtain the NQC parameters (af™,6f™); and finally 
we iterate the procedure until (af™, 6f™) converge at the 
fourth digit. 



C. Determining 05(1^) 

At this stage, the only freedom left in the model is the 
value of fflgC*^)- '^'-'^ explain how we investigated 

the phase difference A^^'-'^^^ (ag(j^)) and used it to de- 
termine ag(i^). Actually, we used a two-pronged approach 
towards studying Acj)^'^^^^. We first studied the Q^{u;) 
function defined by the NR data, and compared it to the 
EOB-predicted one. Then, in a second step, we consid- 
ered the time-domain phase difference Acj)^'^^^^ (t) . 

Let us start by explaining how we used the Qcj{uj) di- 
agnostics to constrain the possible good values of ag(i/). 
Since, as we explained above, we could extract from NR 
data a rather accurate estimate of Q^^{uj), we com- 
pared it to the value (55*^^(aj; a§(i')) predicted, for each 
value of ag(z^), by EOB theory. Such a comparison 
(in the q = 2 case) is illustrated in Fig. 7. The top 
panel of this figure shows the EOB-PN and NR-PN 
differences AQ^ = - Q!j^^^ , where X labels ei- 
ther EOB (for the three indicated values of a^) or NR, 
and Q^^^^ is the 3.5PN-accurate, Taylor-expanded ex- 
pression given by Eq. (60). On the other hand, the 
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FIG. 7. (color online) Using the Qu,{ijj) diagnostics to con- 
strain the good values of 05(1/). The figure refers to the 
case q = 2, v = 2/9. Top panel: difference between ei- 
ther or Q'^^ and the 3.5PN-accurate, Taylor-expanded 
Qti^^^ given by Eq. (60). Bottom panel: the lines show 
the differences AQ^ = Q5°^ - 0^^'^=^ for different val- 
ues of ag. The shaded region exhibits the difference AQ^i — 
gNR,N=5 _ qNR,n=4 ^j^g^g AT = 4, 5 labels two different res- 
olutions, respectively medium and high, of the NR data [34] . 
See text for further details. 



bottom panel displays the EOB— NR difference AQ^ = 
gS°^(w;ag(i/)) - Q^^(w) for five different values of ag. 
In addition, the shaded region represents the NR— NR 
difference AQ„ = Q^^-^=^ - Q^^^''^=^ where where 
N = 5 (respectively N = 4) labels the numerical wave- 
form with the highest (resp. medium) resolution [34]. 
The visual comparisons displayed in Fig. 7 arc made 
quantitative in Table VI, which lists corresponding val- 
ues of the EOB— NR phase difference over the frequency 
interval il/(aji, 0^2) ~ (0.07,0.29) obtained from the inte- 
gral 



TABLE VI. Mass ratio q = 2: phase difference A(j) = 0^"" 
accumulated between frequencies cji = 0.07 and UJ2 
0.29 versus 05 obtained using Eq. (67). 



A(f> [rad] 



-90 
-104 
-106 
-108 
-120 



0.0895 
0.0303 
0.0490 
0.0682 
0.1940 



Ac/) : 



dlnu;(Qr^(a;;a^(j.))-Qf:^^(c.) . (67) 



Note that Muj2 = 0.29 approximately corresponds to the 
merger. These phase differences indicate that a good 
range of values of ag(2/9) is roughly between —104 and 

— 108. Within such a range, Atp remains of the order of 
the NR phasing uncertainty as estimated in Ref. [28, 34] 
by comparing the two resolutions TV = 4 and N = 5. 
Note that the small phase differences corresponding to 

— 108 < ag(2/9) < —104 result from a cancellation be- 
tween positive and negative contributions to the above 
integral. However, a look at Fig 7 shows that within this 
range of Og the nonzero values of AQ^ remain of the or- 
der ±0.05 for most of the integration region. Such range 
of values of AQ^^ is comparable to the numerical uncer- 
tainty on Qi^ (at least) during the inspiral, as illustrated 
by the shaded region in the figure. Note indeed that the 
frequency Mui = 0.1 is reached only 150M before merger 
(cf. bottom left panel of Fig. 9). Note also that the fre- 
quency interval 0.2 < MW22 < 0.3 (where the top panel 
of Fig. 7 shows visible differences) only corresponds to 
the last 25M before merger. [The GW frequency 0.2 ap- 
proximately corresponds to the adiabatic LSO crossing, 
i.e. the end of the quasi-adiabatic inspiral]. 

This analysis based on the Q^j diagnostics selects, for 
each value of the mass ratio {q — 1, 2, 3, 4, 6), a range of 
good values of ag(j/). which then needs to be confirmed 
and refined by directly comparing the time-domain phase 
evolution of the EOB waveform to the NR one. We have 
done such an analysis by considering, for each value of 
ag(i^) within the above range, the phase evolution from 
the beginning of the simulation up to merger, and also 
after merger, during ringdown. The comparison up to 
merger only depends on the choices of t^^j.{v) and ag(i^), 
while the comparison during the subsequent ringdown 
also depends on the choices made in attaching QNM's to 
the NQC-corrected pre-merger signal. The time-domain 
phasing comparison allowed us to close up, for each value 
of V, on a more precisely determined value of ag(i^) (with 
an uncertainty of order unity). Actually, depending on 
the criterion we put on the quality of the EOB/NR phase 
agreement, the resulting best values of ag(i^) are slightly 
different. However, in all the cases we have explored, 
we found that the good, z^-dependent values of ag were 
approximately lying along a straight line. 
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For instance, if we require that the time-domain phase 
difference (after alignment), A(/)^'-'^'^^(t; ag(i/)), remains 
near zero in as flat a manner as possible up to the 
merger, we found that one should choose a,1{v) to (ap- 
proximately) lie on the straight line 

ag"^'(i/) = -123-M43(l-4i/). (68) 

This "flat" ag hnc connects ag^^*((7 = 1) = -123 to 
a%^'^^{q = 6) = —50. However, the price to pay to 
have a very close phase agreement up to merger is that 
the subsequent, somewhat coarse QNM attachment de- 
fined by the current BOB prescriptions, will cause, after 
merger and during ringdown, the EOB-NR phase differ- 
ence A(/)^'-'^^^(t; aQ^^^{v)) to jump to positive values of 
order ~ -f-0.15 rad (more about this below). 

On the other hand, we preferred to look here for an 
"effective" description of the phasing where we allow 
^^EOBNR ioke slightly negative values just before 
merger, but to jump to smaller values ^ -1-0.05 rad after 
merger (see more details below). We can then summa- 
rize our search for such an effective a1{v) by the following 
simple analytical representation of the function a'^{i')'. 

ag(i^) = -122 + 147(1 -41/). (69) 

This is one of the central results of our work, and one of 
the most important pieces in the NR-completion of our 
BOB model. 

In conclusion, we propose to define the NR comple- 
tion of our EOB model by adopting the values (66) 
and (69) for defining, respectively, t^^^{v) and ag(i/). In 
addition, we found that the following QNM-attachment 
choices define a reasonably accurate ringdown comple- 
tion of the EOB waveform: = 5 QNM modes, and 
^comb ^ 0.7M. In the following, we shall illustrate the 
comparison of the EOB multipolar waveform defined by 
these choices to the corresponding NR multipolar wave- 
form. 



D. Values of the NQC parameters (af™, fof") 

Before doing so, let us recall that, for each mass ra- 
tio, we must determine (by iteration) the NQC param- 
eters {ai,bi) defined by the above choices (using given 
NR data). In Table V we list, for the mass ratios 
q (1,2,3,4,6) and for multipoles (2,2), (2,1), (3,3), 
(3,2), the values of the (of™, &f™)'s that define the NQC 
corrections to the bare inspiral-plus-plunge EOB wave- 
form. [When q — 1 there are no entries for £m = (2, 1) 
and (3,3), because these modes are identically zero in 
this case for symmetry reasons.] We will discuss below 
the issue of replacing the information contained in this 
table by i/-dependent fitting formulas. 



E. Effect of the NQC factor on the EOB waveform 

Let us first illustrate how the NQC factor modifies the 
purely inspiral EOB waveform. The q = 1 case is con- 
sidered in Fig. 8: modulus (left panel) and frequency 
(right panel). Similar results are obtained for any other 
mass ratio (see also Ref. [54] for the test- mass limit). 
We show together: (i) the purely inspiral waveform, i.e., 
Eq. (14) without the NQC factor h^^'^'^ (dash-dotted, thin 
line, black online); (ii) the inspiral-|- merger waveform, in- 
cluding the NQC factor (dash-dotted and thick line, blue 
online); (iii) the extended EOB waveform, including the 
ringdown part (thick, solid line, red online); and the NR 
waveform (thin, solid line, black online). As noted al- 
ready in Ref. [54] the most striking feature of this plot 
is that the pure inspiral EOB waveform modulus peaks 
(after alignment as explained in Sec. VF) just ^ 1AM 
before the peak of the NR modulus. On the other hand. 
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FIG. 9. (color online) Comparison between EOB and NR (Zerilli-normalized) waveforms for mass ratios q = 1, 2. Left panels: 
modulus and frequency. Right panels: phase difference (top) and real part (bottom). 



its amplitude is about 20% larger than the NR one®. Note 
that the largish difference in amplitude is very efFectively 
corrected by the NQC factor. In order to reduce the am- 
plitude and displace it to the right we need a NQC factor 
that, near merger, is smaller than one and growing. This 
is what 712 succeeds in doing thanks to its shape, as illus- 
trated in Appendix A. This explains why the values of 
the NQC parameter a^^ are the dominant ones, see Ta- 
ble V. By contrast, if one has to increase the amplitude 
and displace it to the right (as was needed in Ref. [28] 
because of the use of the argument il^/'^ in ^22(2^)), one 
needs a NQC factor which, near merger, is larger than 
one and growing, as, for instance, our ni, Eq. (26a). 



Such a behavior follows from our use oi x = as argument 
in p22ix). As noted in Fig. 2 of Ref. [28], the different choice 
X — f n^/^ (which is however not physically justified during the 
plunge), makes the EOB waveform peak considerably earlier (by 
6.2M) than NR, but with an amplitude much closer to the NR. 
one (~ —0.23% smaller). 



F. Comparison between the I = m = 2 NR and 
EOB waveforms 

Let us now present the results of the comparison be- 
tween the dominant quadrupolar {{i,m) = (2,2)) NR 
waveform, and the corresponding NR-completed EOB 
waveform introduced in this work. For each mass ra- 
tio among q = (1, 2, 3, 4, 6) , Figs. 9-10 compare the EOB 
and NR modulus and frequency (left panels), the real 
parts of the waveforms (right panels, bottom) and the 
phase difference A^^^obnr ^ ^eob „ ^nr (^^.^^^it pan- 
els, top). The vertical dashed line present in all pan- 
els marks the location of the peak of the EOB orbital 
frequency, i^plf^k- These time-domain comparisons are 
done by suitably determining a relative time and phase 
shift between the two phases <l)^^{t^^) and 0f2°^(^^°^)- 
These shifts are estimated by minimizing the time inte- 
gral of the square of the phase difference on a time inter- 
val corresponding to a given frequency interval [wl,^/?]. 
Following Refs. [28, 65] , we perform this waveform align- 
ment on the long inspiral phase. Note that, in doing so, 
we do not enforce the constraint that t^^^ corresponds 



21 




FIG. 10. (color online) Comparison between EOB and NR (Zerilli-normalized) waveforms for mass ratios g = 3, 4, 6. Left 
panels: modulus and frequency. Right panels: phase difference (top) and real part (bottom). 



to i^^pcajj. However, the EOB/NR agreement is so good 
up to merger that such an early-inspiral alignment suc- 
ceeds in realizing, a posteriori, a near coincidence be- 
tween t^^j. and tfjpoak- The right limit of the frequency 
for each mass ratio is u>r = 0.1. The left bounds arc 
cjL = (0.035,0.035,0.035,0.044,0.045). 



These figures indicate an excellent EOB/NR agree- 
ment in phasing and in modulus from the early inspi- 
ral up to merger. The remaining disagreements are well 
within the nominal error bar on numerical data. Ac- 
tually, the only estimate of the numerical error on the 
phasing of these numerical data that is available in the 
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FIG. 11. (color online) Subdominant multipoles, I = 2, m = 1 (left panels) and ^ = m = 3 (right panels) Comparison between 
EOB model and NR (Zerilli-normalized) waveform for mass ratio q = 2. Top: modulus and frequencies. Bottom panels: phase 
differences. 



literature is a rather conservative one that is done by tak- 
ing the difference between the highest and the medium 
resolution. This procedure gives uncertainties that arc 
very small during the inspiral phase (< 0.01 rad) and 
small, though not negligible, in the late plunge phase up 
to merger (^ 0.1 — 0.3 rad, depending on the mass ra- 
tio) [34] . A less conservative NR error estimate might be 
smaller by (at least) a factor two^. Keeping this in mind, 
it is remarkable that our EOB model, with the very sim- 
ple law for a%{v) given in Eq. (69) is able to reproduce 
all numerical data within < 0.06 radians at merger. 

Let us also emphasize the good agreement between the 
moduli around merger, with an excellent global represen- 
tation of its time evolution around merger, and during 
the subsequent ringdown. This is a clear improvement 
with respect to previous works [22, 28, 53] which is due 



® We thank Harald PfcifFer and Luisa Buchman for informing us 
of this more reahstic estimate of the NR errors. 



to a combination of effects coming both from the use of 
an improved analytical EOB model, from a new choice 
of the basis of NQC functions rii, and from the choice 
of an NQC determination point which differs from the 
maximum of the amplitude. Let us also note that, as 
already mentioned, we have, on purpose, chosen effective 
values of ag(z^) causing the phase difference A^^*-*^^^ 
to dip towards negative values ^ —0.05 rad just before 
merger, before jumping towards positive values of order 
+0.05 or +0.1 rad during ringdown. Such a behavior 
ensures a good average phase agreement during the en- 
tire process. Had we instead chosen the slightly different 
"flat" values of ag(j'), Eq. (68), they would have lead to 
a near perfect phase agreement up to merger. However, 
the price for doing so would then have been the presence 
of a larger global phase disagreement (of order ~ +0.15 
rad), due to a positive jump in A</)^°^^^ after merger, 
and during ringdown. We note that such a positive jump 
^ +0.15 rad in A(/)^'^^^'^ is consistent with the study of 
the intrinsic error in Ac/)^*^^^^ coming from the proce- 
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FIG. 12. (color online) Subdominant multipoles, I = 2, m = 1 (left panels) and I = m = <i (right panels) Comparison between 
EOB model and NR (Zerilli-normalized) waveform for mass ratio g = 6. Top: modulus and frequencies. Bottom panels: phase 
differences. 



dure of QNM attachment itself done in Ref. [28]. This 
indicates that more work should be devoted towards im- 
proving the current EOB technique for attaching QNM's 
onto the inspiral-plus-plunge waveform. 



G. Subdominant multipoles 

Up to now, our study has only considered the dom- 
inant quadrupolar I — m = 2 waveform. Let us now 
compare some of the subdominant multipolar waveforms. 
We consider here the i — 2. m = 1 and £ = m = 3 sub- 
dominant waveforms, for the two mass ratios q = 2 and 
g = 6 (similar results were obtained for q = 3 and (7 = 4). 
We limit ourselves to such a partial comparison here to 
show the capability of the EOB model, as it was defined 
above, to get the main characteristics of the subdominant 
multipoles, without introducing ad hoc modifications, or 
tuning further parameters. At the end of this section 
we will also mention some results for the £ = 3, m = 2 



multipole. 

In Figs. 11-12 we compare, for the two mass ratios 
q = 2 and q = 6, the NR and EOB frequency and mod- 
ulus for the two subdominant multipoles £ = 2, m — 1 
and ^ = m = 3. We use the same matching interval 
as for the £ = m = 2 mode, i.e. A™^*'^^ = 0.7M, and 
the same number of QNM modes, i.e. N = 5. Note the 
good agreement of the moduli in all cases, both up to 
merger, and during ringdown [In the A21, q = 2 case the 
multiple crossings between the NR and EOB moduli may 
be due to inaccuracies in the NR waveform.] Note also 
the good agreement, up to merger, of the frequencies, in 
all cases, and the good agreement of the frequency of the 
(3, 3) mode after merger, and during ringdown. The only 
case which is slightly less successful is the discrepancy be- 
tween the EOB frequency and the NR frequency in the 
£ = 2, m = I case for both mass ratios (compare with 
Ref. [28], but note we have not introduced here any ad 
hoc treatment of the the £ ~ 2, m = 1 case.) Namely, the 
EOB frequency of the (2, 1) mode shoots up, just after 
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subdominant multipoles are determined to be 
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FIG. 13. (color online) Subdominant multipole ^ = 3, m = 2: 
phase difference for q = 1 (top panel) and g = 6 (bottom 
panel). The vertical dashed lines mark the f^pfak crossing 
time. 



merger , a bit faster than its NR counterpart. In turn, 
such a frequency difference builds up a phase difference 
after merger. This is illustrated in the bottom panels 
of the figure, which show the phase difference A^ff*^'^^ 
(left) and A^g^^^^ (right) as a function of time dur- 
ing the entire simulation. Note that the dephasing is 
remarkably small up to merger for both multipoles, and 
then accumulates a dephasing A0ff ^'^^ ^0.5 rad (and 
Ac^g^^'^^ 0.15 rad) during the ringdown. 

Let us emphasize that the phase difference 
A0f °^NR(^) plotted in the bottom panels of Figs. 11-12 
has been computed without introducing any new arbi- 
trariness, neither in time, nor in phase, in comparing 
the two phase evolutions. Indeed, the least-squares 
alignment procedure of the NR and BOB dominant 
(2,2) waveforms has determined both a shift in time, 
say r22, and a phase shift, say Q!22, connecting them. 
The time shift T22 determines the (a priori unknown) 
connection between the two time variables t^^ and 
i^*-*^, and should therefore be used in comparing the 
time evolutions of all the other physical quantities, and 
in particular the subdominant multipoles. The case of 
the phase shift a22 is similar, but with a difference. 
Indeed, in our case (with a common, preferred z axis 
given by the total angular momentum of the sytem) 
the only a priori unknown angular difference between 
NR and EOB is a rotational shift, by some angle /3, 
connecting the NR basis of tensorial spherical harmonics 
to the corresponding EOB basis. This common angle 
P then introduces a phase shift in all the various ^m 
multipoles simply given by 



Oitm = m/3, 



(70) 



independently of £. As this result applies in particular to 
a22 (which is determined modulo 27r by the alignment of 
the (2, 2) waveforms), we sec that the phase shifts in the 



■ a22 modulo rrnr 



(71) 



In addition to this phase shift, there might be extra phase 
shifts due to the use of different conventions in defin- 
ing the phase of the tensorial spherical harmonics. Such 
phase conventions differ at most by multiples of 7r/2, cor- 
responding to powers of i. In other words, we can always 
write that a^m = y "^22 modulo 7r/2, which is sufficient 
for unambiguously computing A(t)f^^^^ for all subdom- 
inant multipoles. This absence of phase-shift ambiguity 
in A<j>f^^^^ makes it all the more remarkable that, in 
the (2, 1) case, the phase difference At/iff ^^'^ plotted in 
Fig. 11 (for q = 2) and Fig. 12 (for q = 6) stays very 
small up to merger. 

Let us finally comment on Fig. 13, were we show the 
phase difference one gets for the £ = 3, m = 2 multi- 
pole, for the two representative cases q = I (top panel) 
and q = 6 (bottom panel). The figure, again, illus- 
trates a rather good consistency between EOB and NR 
up to merger. The differences after merger are mostly 
due to our simplified description of the ringdown (see 
Appendix A of Ref. [28] for a detailed analysis of the 
structure of the (3, 2) ringdown waveform). 

We leave to future work a more detailed analysis of the 
subdominant multipoles, and the investigation of pos- 
sible ways of improving their EOB representation, in 
case the slight dephasing exhibited in Figs. 11-12 for the 
{£, m) = (2, 1) multipole happens to significantly degrade 
the faithfulness of the complete EOB waveform (summed 
over all multipoles). 



VI. STRUCTURE OF THE EOBNR RADIAL 
POTENTIAL A{u) AND ITS CONNECTION 
WITH OTHER RESULTS 



One of the most important nonperturbative dynami- 
cal knowledge acquired in this work by comparing EOB 
predictions to the Caltech-Cornell-CITA simulations con- 
cerns the function A(u; v). We recall that A{u] v) is the 
main radial potential of the EOB Hamiltonian, and repre- 
sents the time-time component of the effective EOB met- 
ric: A{u]v) = — (7§o(i?). In the test-mass limit, — )■ 0, 
the effective metric is the Schwarzschild metric, so that 
\mY^^oA{u;v) = 1 - 2u = 1 - 2GM/{Rc^). We saw 
above that NR data selected, in the strong-field domain, 
an A function given by Eq. (10) with Og = 23.5 and 
a1{v) = —122 -I- 147(1 — Av). Let us now discuss some 
properties of this NR-informed EOB potential (or simply 
EOBNR potential) and its connection with other relevant 
results. 
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FIG. 14. (color online) Contrasting various estimates of the 
A{u; v) function, in the equal-mass case, v = 0.25. The plot 
shows the IPN, 2PN and 3PN Taylor-expanded versions of 
A{u; 0.25); its 3PN-accurate Pade resummed form as well as 
the EOBNR one (5PN accurate with logarithmic terms and 
al = 23.5, ai = -122). 



A. Global shape of A^'^^^'^{u; v) as a function of u, 
and comparison virith previous purely analytical 
estimates 

As a first orientation, we contrast in Fig. 14 vari- 
ous estimates of the function A{u] v) in the equal-mass 
case, i.e. v = 0.25. Our NR-informed estimate (5PN- 
log-Pade resummed and with a§ = 23.5 and a%[v) = 
— 122 -t- 147(1 — Av)) is shown as a thick solid line (red 
online), i.e. the second line from the top. The dashed 
bottom line represents the IPN-accurate estimate of A, 
which happens to coincide with the simple Schwarzschild- 
metric result A^*'^(m) = 1 - 2u. [Indeed, in Eq. (7) 
there are no terms of order w^, corresponding to the IPN 
level.] The thicker dashed line just above this IPN esti- 
mate represents the Taylor-expanded 2PN estimate, i.e. 
Eq. (7) taken up to the term 0{u^) included. The up- 
per dashed line represents the Taylor-expanded 3PN es- 
timate of A(u\v), as given by Eq. (7) up to the term 
Oiu^) included. Finally, the thin solid line (black on- 
line) just below the NR-completed 5PNlog Pade curve is 
the Pade-resummed estimate of the analytically known 
3PN result, which was proposed by Damour, Jaranowski 
and Schafer [3] in 2000, i.e. five years before NR simula- 
tions started yielding information about the strong-field 
dynamics of binary black holes. It is remarkable that the 
latter simple 3PN-Padc estimate is rather close to the 
best current NR-informed estimate: (i) it is numerically 
quite close to it if one considers values m < 0.3 which 
are already beyond the last stable orbit, and therefore 
are crossed during the plunge; and (ii) even in the very 
strong field domain 0.3 < u < 0.6 (where the merger 
occurs) the 3PN-Pade estimate is a much better approx- 
imation to A^^^^^{u]v) than any of its standard PN 



approximants. This closeness explains the success of the 
simple Padeed 3PN A function in agreeing with several 
recent NR studies of dynamical aspects of close black hole 
binaries [26, 33], and confirms the effectiveness of using 
Pade approximants to improve the strong-field behavior 
of Taylor approximants. 



B. Detailed study of the //-dependence of 

^^OSNi^(u;i.) 

The comparison of the previous subsection has indi- 
cated that an accurate description of the gravitational 
wave emission of coalescing binary black holes requires 
a very precise determination of the shape of A{u] v) in 
the very strong-field domain u > 0.3 (i.e. R < 3GM/c^). 
Let us zoom on the detailed shape of the A function in 
the strong-field domain by focusing on the properties of 
the associated a function, defined by writing 

A{u;iy) = l-2u + i^a{u;i^). (72) 

The Taylor expansion of this small-a function starts as 

'94 41 



a{u; v) = 2u'^ + 
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Ti^ + 0{u''\nu). (73) 



Note that the v dependence of a{u; v) is only contained 
in the 0{u^hiu) remainder term. In order to zoom on 
the V dependence of a(u; v) it is then useful, following 
Ref. [47], to normalize the a function by its LO PN be- 



havior, a^^^{u]v) 
function defined as 



= 2u'^, i.e. to consider the a{u;v) 



a{u; v) 



a{u\ v) _ A(u\ v) — {1 — 2u) 



2u3 



2vu^ 



(74) 



In the upper panel of Fig. 15, we plot the values of the 
EOBNR a{u\ v) functions for the values of v correspond- 
ing to the five mass ratios we used in our EOB /NR com- 
parisons above, namely q = 1,2,3,4,6, as well as the 
EOBNR predicted d curves corresponding to q ~ 10, 
to q = 100 and also io q = oo, i.e. to the v = 
q/{q + 1)2 ^ limit of a^'^^^^ {u; u) . The (red on- 
line) round markers on the curves indicate the EOB- 
defined, light-ring locations, i.e. the solutions of the 
equation {v?A{u))' = (see Table VII for the precise 
numbers). In addition, we have also indicated the re- 
cently derived (GSF-computed) "exact" value of the limit 
lim^_j.o a(u; v) [47] (using their best analytical fit). In the 
bottom panel of Fig. 15 we plot the corresponding val- 
ues of the products i/aEOBNR(j^. ^) ^ 2vu^a^^^^^{u] u), 
i.e. the corresponding differences of ^^'-'^^^(u; v) away 
from its test-mass limit, i.e. A^'^^^^ [u; v) - A^^'^'^lu), 
where A^'^'^^iu) = I - 2u = lim^_,o ^^^^^^(u; v). This 
shows again how the physics of the GW emission by 
coalescing black hole binaries depends on fine features 
in the A potential. Note how, as v decreases, a(u; v) 
monotononically increases, in a way which is qualitatively 
compatible with the shape of the limiting GSF result 
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FIG. 15. (color online) Top panel: Behavior of the EOBNR 
a{u) function defined in Eq. (74) with = 23.5 and 05 — 
-122 + 147(1 - 4:1/). The red line shows the 1/ = function 
as obtained from the fit of GSF data [47]. Bottom panel: the 
difference A^°'^n^(m; u) - A^'^'^'^iu) with A^'''^'~{u) = 1 - 2u. 
For each value of u, the marker indicates the EOB-defined 
adiabatic hght ring location. 



a(u; 0) = lun^^o o-i^'i i^)- [The latter limiting GSF shape 
has a singularity at u = 1/3, which is probably smoothed 
out by higher-order corrections in v around ly ~ 0. See 
[47] for a detailed discussion of the origin of this singu- 
larity, and its probable fictitious character.] Though the 
J/ — > limit of a^*^^^^ (u; ly) (which is a polynomial in u, 
with logarithmic coefficients) does not coincide with the 
exact 0{v) GSF result, it stays quite close to it up to 
u < 0.2. Note also, on the bottom panel, how the behav- 
ior of the corresponding contribution to the A potential, 
i.e. the product i'a{u\i/)^ seems to tend continuously 
(though maybe not uniformly) towards zero as u Q. 
This bottom panel suggests that the g = 10 case should 
be thought of as belonging to the class of the normal 
comparable-mass cases q = 0{1). One needs g's of order 
at least 0(100) to belong to the class of extreme-mass- 
ratio binaries. The EOBNR potential derived here has 
anyway been tuned to the physics of comparable-mass 
binaries with 1 < g < 6. As we knew (from Rcf. [47]) 
that the v ^ Q limit of the (exact) A potential was (prob- 
ably) mildly singular, and as we are mainly interested in 
describing the physics of comparable-mass systems, we 
did not attempt to incorporate in the A function too 
much of the information contained in its — >■ 0, GSF 
limit. In our work above, we only incorporated some in- 
formation about the — >■ limit of the 4PN coefficient 
lim^^o fl5(i^)- But, as we shall discuss next, this was 
mainly done as a practical way of reducing the number 
of unknowns to be fitted to NR data. 



TABLE VII. EOB-defined light-ring (LR) and last-stable- 
orbit (LSO) locations for al = 23.5 and al = -122 + 147(1 - 
4v). 
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0.25 


1.8245 


0.5481 


4.4904 


0.2227 


2 


0.2 


1.9585 


0.5106 


4.6847 


0.2135 


3 


0.1875 


2.133 


0.4688 


4.9210 


0.2032 


4 


0.1600 


2.2665 


0.4412 


5.0969 


0.1962 


6 


0.1224 


2.442 


0.4095 


5.3238 


0.1878 


10 


0.0826 


2.633 


0.3812 


5.5519 


0.1801 


00 





3.0000 


0.3 


6.000 


0.16 



C. On the "equivalence classes" of the A(u) 
potential 

References [6, 22] found, for the q=\ case, that there 
was a strong degeneracy between the two parameters en- 
tering a 5PN-accurate Fade representation of the A func- 
tion, say a\ and Og. This was confirmed for other values 
of q in Ref. [28]. This finding leads to the idea that the 
good values of 05 and a% can be organized in "equiva- 
lence classes" of quasi-interchangeable values of the pairs 
(05, a§). An explicit way of constructing these equiva- 
lence classes was indicated in [35]: it consists in defin- 
ing the equivalence class of some given pair (a^ , Og ) 
as the set of pairs (a§,ag) such that the u-derivative 
A'(u; v\ a§, flg) of the A function, evaluated at some fidu- 
cial strong-field point, say Wb (the value ~ 0.215 was 
suggested there), takes the same value at (og, Og) and at 

{a\ , flg ) . In equations 

A'(ub;z.;a^,ae)=A'(^.b;^^;a^^°\a^^"^) , (75) 

or, equivalently, 

a!{ub\v-,al,a%) = a!{ub\v\a%^^'' ,a%^^^) . (76) 

When working, as we do here, with the normalized func- 
tion a(u;j^), we could alternatively define these equiva- 
lence classes as level sets (in the space of pairs (a§,a§)) 
of a'(u6; i^; 05, a§), or even, simply, of a(ub; z^; 05, a§). 
Evidently, all those possible "definitions" lead (when 
one changes the fiducial value Uh, and/or the consid- 
ered function a', a', a, etc.) to different equivalence 
classes. However, because of the properties of the A 
function, one checks that, as long as one bases one's 
definition on the value of A or some related function 
in the strong-field region, this leads, to a good approxi- 
mation, to a numerically rather well-defined equivalence 
class of (a§,ag) pairs. This is illustrated in Fig. 16. 
This figure shows (for the case q = 1) that our NR- 

tuned preferred values {al^^\ flg*'*^'') — (23.5, —122) define 
a i'a{u) function which can be very nearly reproduced 
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by using other pairs of (ogjOg) values, namely (0,220), 
(5,125), or (10,40). The upper panel shows together, 
versus u, the functions i^a(ti; Og, Og), for g = 1, i.e. 
V = 0.25, and for the four different pairs of parameters 
values (a§,ag) = (23.5, -122), (10, 40), (5, 125), (0, 220). 
The plot illustrates that these five different functions 
are indishtinguishable by eye. The bottom panel of 
the figure zooms on the differences away from our stan- 
dard choice (a5^°\a6^°^) = (23.5,-122), i.e. it plots 
;/Aa(u; a§, Og) = A(u; z/; a§, ag) — A(u; i^; 23.5, —122). 
For any choice of the parameters, these differences are of 
the order lO"**. Note that we have not used, here, any 
precise, level-set type, criterion for selecting the pairs 
equivalent to our preferred value, but we have selected 
them by simple trial and error, until we could reduce the 
(maximum) difference to the smallest level we could find . 
This smallest level was O(10~'*). The reason why such a 
level of deviation is small enough for our purpose can be 
seen by turning back to our analysis above, when we were 
fixing the fiducial value a§ = 23.5, and then tuning the 
value of a% for the EOB phasing to best agree with the NR 
one. In that case (as is clear from the number of digits we 
were giving above for Og) we found that the "good" val- 
ues of a§ were determined, roughly, within an uncertainty 
8a% = 0(1). Such an uncertainty on the good value of a% 
(for the fixed 05 = 23.5) entails a corresponding uncer- 
tainty on the value of the function ^(u;a§,ag) of order 
5A{u;a's^,a%) ^ dA{u;a's^,a%)/da%. The latter quantity is 
found to increase with u, and to reach a value of order 
0.8 X lO^'' when u takes the light-ring value ulr — 0.55 
(for q = 1). In conclusion, a possible variation in the 
A{u) function of Loo norm ~ 10~^, for < u < ulr, 
is a reasonable way of defining the equivalence class of 
A{u), and Fig. 16 shows that one can indeed, starting 

from our values {a'i}'^\ Og*""'') = (23.5, —122), find a (rela- 
tively thin) strip of values of (ag, a§) along which the 5PN 
Faded function A^*-*^^^ (u; a§, ag; v) stays within such an 
equivalence class. 

Though here wc focus only on the q — \ case, similar 
classes of equivalence of a functions exist for any mass ra- 
tio. In summary, this exercize confirms that we were jus- 
tified in a priori fixing the value of 05 . Finally, the impor- 
tant fact is that NR data allow one to directly determine 
the A[u] v) function itself, essentially independently of 
the chosen "representative" (a§,a§) within some equiv- 
alence strip in the (a5,ag) plane. This determination 
of the A^'^^^^ {u\ v) function is exemplified on Fig. 15 
(keeping in mind the invisible deviations plotted in the 
upper panel of Fig. 16). 



D. Comparison between the present determination 
oi A^'^™'^{u-v) (5PN with logs), with previous 
estimates (5PN without logs) due to Damour and 
Nagar (2009) and Pan et al. (2011) 

The present work is the first EOB work to include log- 
arithmic terms in a comparison with NR data. Let us 
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FIG. 16. (color online) Elements of the equivalence class of 
£(05, ag) functions for g = 1. The bottom panel shows the 
fractional difference with our favourite choice a% = 23.5, ag = 
-122. 



now compare our final NR-aided determination of such a 
A function (with logarithmic terms) to the 5PN-accurate 
A function [without logarithmic terms) used in previous 
EOB works [6, 22, 27, 28]. In particular, Ref. [22], using 
a 5PN-accurate A{a^,ag; v) function (without logs), ex- 
ploited a previous version of the q ^ 1 Caltech-Cornell- 
CITA numerical waveform to find a banana-like region 
of good values in the [a'^^af) plane such that the phase 
difference between EOB and NR waveform through in- 
spiral, plunge and merger was < 0.02 rad. The values 
a% = —6.37 and a% = 50 lied in the middle of this 
good region and have been used extensively in subse- 
quent EOB work [26, 27, 31, 42]. [By contrast Ref. [22] 
actually used the values a§ = and a§ = — 20 which lie 
on the boundary of the good region.] Then, the analog, 
banana-shaped equivalence classes in the (a5,a§) plane 
corresponding to other values of q were first investigated 
in Ref. [28]. [The latter reference basically used the same 
conceptual structure as Ref. [22] with some technical dif- 
ferences.] Ref. [28] found a very good agreement between 
EOB and NR waveforms with an A function defined by 
the following choices 

az{u) = -5.828 - 143.5;/ + 447^^ (77) 
ag = 184 . (78) 

In Fig. 17 we consider the two mass ratios g = 1 and g = 6 
and for each mass ratio we compare three different a{u) 
curves, namely: (i) the log-containing 5PN-accurate one 
determined in this work ("EOBNRlog" with ag = -)-23.5, 
and al = -122 + 147(1 - Av)); (ii) the log-less 5PN- 
accurate one of [22] (with a§ = —6.3 and a% = 50); and, 
(in) the log-less 5PN-accurate one of [28], see Eq. (77). 
The figure shows that while the three different analyti- 
cal descriptions seem to be visually close for the equal- 
mass case, g = 1, they exhibit visible differences in the 
(7 = 6 case. However, we have seen above that only differ- 
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FIG. 17. (color online) Comparing a functions for different, 
5PN-accurate, EOBNR-completed models. The markers in- 
dicate the location of the EOB-defined adiabatic light-ring for 
each curve. See text for explanations. 



ences of order 10 ^ in the A function can be considered 
as being negligibly small. When computing the differ- 
ences AA^{u; u) = A^~ ^EOBNRiog ^j. ^j^g j^^^^gig 

X = DN2009, Ref. [22] and X = Pan et al, Ref. [28], one 
finds that, for q = 1, AA-^ (u) is a monotonically decreas- 
ing function of u which reaches values of order ~ —0.004 
for X ^ DN2009 and ~ -0.0025 for X = Pan et al when 
u ~ 0.5, i.e., close to the corresponding adiabatic light- 
ring position. Such differences are therefore quite signif- 
icant on the 10~^ scale of our present equivalence classes 
of A functions. In the q = 6 case the corresponding dif- 
ferences taken at m ~ 0.4, close to the adiabatic light-ring 
position, are ~ -0.006 for X = DN2009 and ~ -h0.003 
for X = Pan et al.. Again these differences are quite 
significant. Note however that for u < 0.3 the log-less 
model of [28], Eq. (77), (which had been tuned to the 
same q ~ Q NR data as ours) stays quite close to our 
present log-containing model {A A = 2 x lO^"'). The 
present analysis of these differences in the A function 
shows that the introduction of logarithmic contributions 
in the A function cannot be reabsorbed by tuning other 
structures of the EOB model. As we know, from ana- 
lytical PN work, that these logarithmic contributions do 
exist, we conclude that it is necessary to include them, 
and therefore to abandon previous log-less versions of the 
EOB Hamiltonian in favour of the type of improved EOB 
model presented in this work. 



VII. EXTENSION OF THE MODEL BY 
ANALYTIC CONTINUATION IN v 

In the present work, we have used a discrete sample 
of numerical simulations to complete an EOB model, no- 



TABLE VIII. Fits of Zerilli-normaUzed multipolar quanti- 
ties (amphtude, frequencies and derivatives) extracted at the 
iejrtr(j^) as function of v. Each quantity is fitted to a quadratic 
polynomial of the form femi^) = dTv'^ -f c'f'v -f 4'". For 
the amplitude and its derivatives the full leading-order in- 
dependence {yct,A^t(y), see Eq. 20) is factorized before fitting. 
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tably through the use of suitable, NR-fitted NQC correc- 
tions. In order to be able to compute the predictions of 
such a NR-completed EOB model for arbitrary values of 
I/, we need to fix a procedure for computing the six NQC 
parameters, {af"^ (ly) , bj"^ {v)) , as continuous functions of 
i/. [The remaining defining elements of the EOB model, 
notably and Og were already given as functions of v.] 
One can think of two different ways of continuously ex- 
tending the definition of the present EOB model to any 
value of v: first, one can interpolate the discrete sam- 
ple of (a^™, bf^) values of the NQC parameters that we 
obtained (from the five numerical simulations with q = 
1, 2, 3, 4, 6) by fitting them to, say, quadratic polynomials 
in v; second, one can instead fit the original NR-extracted 
numerical values of {Ai„i, Ae,n, Aim,u}im,u}im:U)em) to 
quadratic polynomials in and then, for any given v, 
determine {af"^ {v) , bf"^ (v)) with the iterative procedure 
described above. We have explored in detail both pro- 
cedures. The first one, i.e fitting the end parameters 
(af™,&f'") needed to compute an EOB waveform would 
clearly be a faster way to compute, for any v, a corre- 
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TABLE IX. Fits of the ^ = 2 NQC parameters (af 
function of v. Each quantity is fitted to a quadratic polyno- 
mial of the form fi,n (v) 



C2 V + Ci U + Co . 
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spending EOB waveform. Indeed, this approach does not 
require any iteration procedure. We found that the fit- 
ted (af 6f'")(^)'s give very accurate results for the 
£ = m = 2 waveform, and sufficiently accurate results 
for the £ = 2, m = 1 one. Unfortunately, our simple 
quadratic fits did not perform as well for the other mul- 
tipoles. Therefore if one is mainly interested in having 
accurate models for the £ = 2 multipoles (and is ready to 
have less accurate models for the other values of £), the 
simple quadratic fits of the (af"*, bj"^) values listed in Ta- 
ble IX can be used to efficiently generate EOB waveforms 
for any v. 

On the other hand, if one wishes to be sure to re- 
produce, for any multipole, and any z/, EOB waveforms 
which have the same quality as the ones we showed in 
the figures above, one should either construct better fits 
for the NQC parameters (checking their compatibility 
with fits of the extracted NR data), or, instead, start di- 
rectly from the quadratic fits to the originally extracted 
NR data that we give in Table VIII. Contrarily to the 
fits of the (fli, bi) mentioned above (which relied only on 
the q = 1,2,3,4,6 data), we have done quadratic fits of 
{Aim, Ae^,uJera,i^evi,'-^eni) to six numerical results, 
namely the Caltech- Cornell- CITA q = (1,2,3,4,6) data 
together with the q = oo data of Ref . [42] . 

Given these fits, one then needs to solve for the NQC 
parameters. Actually, such a procedure is simplified by 
the fact that, as we said, the quadratic fits for the ap(;/)'s 
(which are the only NQC parameters which need to be 
re- inserted in the flux) can be used from the start, so 
that, contrary to the general case, one can get the needed 
values of the other NQC parameters in one go, without 
having to iterate the procedure. 

Figure 18 illustrates the performances of the two dif- 
ferent fitting procedures. The figure refers to mass ratio 
q ~ 2 only (equivalent results are found for the other 
mass ratios) and shows the following triple comparison 
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FIG. 18. (color online) Testing two possible fitting strategies 
to continuously extend the discrete sample of NQC parame- 



modulus and frequencies. See text for explanations 



ters (a'"', 6,™) to any value of v. comparison between £ — 2 



for £ = m = 2 (top) and ^ = 2. m 1 (bottom) between: 
(i) the NR waveform frequency and modulus; (ii) the 
EOB waveform frequency and modulus obtained using 
the fits (af™(i/), 6f™(i^); (in) EOB waveform frequency 
and modulus obtained by fitting the numerical data ex- 
tracted at t^^^, determining the NQC parameters in the 
usual way and then iterating a few times. 



VIII. CONCLUSIONS 

We have improved the EOB description of nonspinning 
coalescing black hole binaries by incorporating several 
recent analytical advances, namely: 

(i) 4PN and 5PN logarithmic contributions to the con- 
servative dynamics [35-38]; 
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(ii) the 0{iy) 4PN nonlogarithmic contribution to the 
conservative dynamics [36,38-40]; 

(iii) rcsummed horizon-absorption contributions to an- 
gular momentum loss [41, 42]; 

(iv) the radial component of the radiation reaction force 
implied by consistency with the azimuthal one [43] ; 

(v) an additional 3.5PN contribution to the phase of of 
the (factorized [5, 8, 9]) quadrupolar waveform [44]. 

Moreover, we have introduced new features in the EOB 
formalism, namely: 

(a) a Pade resummation of the additional tail phases 
Sim of the factorized EOB waveform; 

(b) a new way of matching the EOB waveform to the 
NR one by mapping the EOB time when the orbital 
frequency reaches a maximum ifjpoai? a specifi- 
cally chosen (z/-dependent) NR time t^^j.{i') around 
merger, Eq. (53). More precisely, we impose [by 
using six next-to-quasi-circular (NQC) parameters] 
a contact between the amplitudes and the fre- 
quencies of the NR and EOB waveforms at the NR 
instant t^^^{i/) which corresponds to top^ak- 

We have extracted new information from the NR data, 
namely 

(c) We showed how to extract from NR (curvature) 
phasing data the function Q^^{uj) = uj^/uj which 
is an intrinsic measure of the phase evolution. We 
have given an explicit representation of the func- 
tion Q^^{uj] q), for q = (1,2,3,4,6), in terms of 
some fitting coefficients (see Eqs. (59), (62) and Ta- 
ble III). 

(d) We extracted data on the NR amplitude and fre- 
quency, together with their first two derivatives, at 
the specific I'-dependent NR time t^^^^i'), which 
is located a little bit after the maximum of the 
quadrupolar waveform amplitude. We gave fitting 
formulas for the i^-dependcncc of those quantities 
for several multipoles, see Table VIII. 

Using such nonperturbative information from NR data, 
we showed how to complete the EOB model by: 

1. Constraining the value of the main EOB radial po- 
tential, i.e. the A(u; v) function; and 

2. Determining the coefficients entering the NQC cor- 
rection factor Eq. (25). 

Among these results, we think that the new expression 
of the NR-tuned A function, containing logarithms, is 
more refined and more accurate than its previous deter- 
minations [22, 28]. Let us recall that, as in previous 
work, the A function is parametrized in terms of co- 
efficients, here called (ag,ag), entering a certain Pade 



approximant,AP'^'i<'(u; i/; a§,ag), Eq. (10). Then NR 
data were used to constrain these parameters. We have 
delineated the reason why the two parameters (a|,a§) 
entering the Pade definition of A^'^'^®(-u; i/; a^.a^) are de- 
generate by giving a definition of equivalence classes of 
the pairs (a§, Og) in terms of some L^o norm of the A{u) 
function. We have determined a good NR-tuned A func- 
tion by assuming a fixed value of (ag = 23.5 as sug- 
gested by recent GSF results [38, 47]), and by then tuning 
the remaining parameter a^^u). We found that a'^ii') can 
be simply represented as a linear function of v 

ag(i^) = -122 + 147(1 -4i/) . (79) 

We think that the resulting function of u and 
^EOBNR(y. ^) ^ AP^dc^y. ^. 23.5,-122 + 147(1 - 4i/)), 
yields an accurate representation of the A{u] v) function 
itself, independently of the way it was obtained. More- 
over, we find remarkable that the good value of A{u; v) 
could be obtained already by considering only the inspiral 
phasing (before the LSO crossing) and was then checked 
to yield (together with the NR-determined NQC correc- 
tions) an excellent phasing agreement up to merger. 

We have presented our improved EOB model in a self 
contained manner so as to allow interested readers to 
generate for themselves all our EOB results. We in- 
tend to make available soon a public version of our EOB 
codes. In view of the new physics that we have in- 
cluded in our EOB model, and of its excellent perfor- 
mance (obtained without introducing any ad hoc param- 
eters) against the very accurate Caltech- Cornell- CITA 
data, we recommend to use this new EOB model in fu- 
ture EOB works (in particular in extensions to spinning 
and/or tidally interacting systems). 
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Appendix A: On the computation of r 

In the definition of the NQC correction 712, Eq. (26b), 
we used for the the second time derivative of the relative 
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FIG. 19. (color online) Top: computation of r with finite 
differencing and analytical approximations and comparison 
with (r)'''^ Bottom: effect on the NQC basis vector n2. The 
figure refers to g = 1 with the choices 05 = 23.5 and Og = 
— 122. See text for discussion. 



separation r the quantity {r)^^\ which is the value of 
r along the conservative dynamics, i.e. neglecting the 
contributions proportional to This choice is made for 
efficiency's sake, since (fy^ is easy to compute along the 
dynamics. In spite of the neglect of J-^ in its computation, 
(f)(0) does represent an allowed NQC correction because 
it vanishes (together with r and the exact value of f) in 
the circular limit (see below). 

For completeness, let us discuss here how to compute 
a more exact value of f along the dynamics and how the 
result differs from (f ) '^^^ . Let us first recall that along the 
EOB equations of motion r is, at any moment, a function 
of the phase space variables: r — r{r{t),p^{t),pr^{t)). 
Therefore, its total time derivative is the sum of three 
partial contributions 



dr 
dr 



dr 
dpr. 



Pr, 



dr 

dp, 



-Pv 



(Al) 



Using the other EOB equations of motion, this equation 
reads explicitly 



dr 
dr 



dr 
dpr. 



dH^ 



EOB 



9?" 



dr p 
dpv 



(A2) 



where dU^o^jdr, = {Aj Bfl-^dH^o^j dr . 

By definition, the circular dynamics limit corresponds 
to setting r = = p^. and dH-^o^jdr = 0. One then 
sees that, along the circular dynamics, one has also Jv^ oc 
Pr. = 0, and (using r = C{r,pr,,p^)pr,) dr/dp^ cc Pr, = 
0. As a consequence, both r and (f )^°^ , defined by setting 
to zero the contributions proportional to i.e. 



(f)(0) ^ ^% dHEOB ^ 

dr dpr, dr^ 



(A3) 



TABLE X. Fits of Zerilli-normalized multipolar quantities 
of the numerical waveforms (modulus, frequency and their 
derivatives) measured at the peak of each multipole ("max- 
ima") as function of v. Each quantity is fitted to a quadratic 
polynomial of the form /f™(j/) = ci""!^^ + ci'^u + c^™. For the 
modulus and its derivatives the leading order //-dependence 
(Eq. 20) is factorized before fitting. 
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vanish in the circular dynamics approximation. This 
shows that wc can use cither the exact r or its "geodesic" 
approximation (f)'^°) to define the second element of the 
"NQC basis", 712 = r/{rn^). 

When using the definition na = (f)(°V(^^^). Eq. (A3) 
allows one to compute immediately n2 along the dynam- 
ics. By contrast, if one wished to use the definition 
n'2 = f/{r^l'^), a comphcation arises. Indeed, as contri- 
butions proportional to J>, and appear on the r.h.s. 
of Eq. (A2), and as these contain the squared modu- 
lus of the NQC factor (i.e., for each multipole, a factor 
1 1 + of-j^nj p) we see that n^ocr now appears on both 
sides of Eq. (A2). 

Schematically, defining ^ 
the structure 



r,Pr,,Pv}- Eq. (A2) has 

r = a(C) + 6(0-Fr. it r) + c(e) r) (A4) 

which only gives and implicit equation for determining 
the exact r along the dynamics. We can however get 
an explicit expression for f by an iterative procedure. 
Inserting f^*'-' as lowest order approximation on the r.h.s. 
of Eq. (A4) defines an improved value, say (f)^^) for f, 
namely 



f:(l) 



+ 6(e)-^r..(e,rW) + c(OJ-^(^,fW). (A5) 
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NQC corrections to the energy flux on the phasing for q — 6. 



FIG. 20. (color online) Mass ratio q — 6: EOB waveform 
(frequency and modulus) obtained by determining NQC cor- 
rections from NR data extracted at t^^^ ^^^y. instead of t^^^. 



By iterating the procedure once more, we then get 

= a(C) + 6(0-^..(f,r-(i)) + c(C)J-^(^,f(i)). (A6) 

The result (A6) leads to a sufBciently accurate computa- 
tion of f up to merger, as illustrated in the top panel of 
Fig. 19. However, the recursive presence of the flux in this 
iteration substantially increases (by approximately a fac- 
tor 4) the computational time needed to produce an EOB 
waveform. This is why we prefer to use rt2 = f*^°-'/(rr2^). 
Anyway, as Fig. 19 shows, n2 and n2 are numerically 
quite similar. In view of the arguments above their dif- 
ferences are essentially absorbed in a redefinition of the 
coefficients a^. 



Appendix B: NQC factor determined using NR data 



at t 



NR 

A22 peak 



In the text, we argued that it was advantageous to 
determine NQC corrections by matching the EOB wave- 
form (considered at p|?a]j) to the NR waveform consid- 
ered at the time t^^j.. Let us illustrate here (see Fig. 20) 
in the case q ~ 6 the slightly different (but significantly 
worsened) EOB waveform obtained when one instead 
matches the £ = m = 2 EOB waveform (considered at 
time ^npcak) ^° ^^e NR waveform considered at the time 
peak (^^ done in all previous EOB works). Fig- 
ure 20 uses as before 6 NQC corrections and the value 
ag = -122+147(1 - 4 X 6/49) = -47. However, the NR 
extraction point, which is also used as NQC determina- 
tion point, is now t^f^ p^^i,. 

The fits of the vector of NR quantities 
iAZ^^Z,u;Z,d^fr^.^Z) now measured at the 



location of the maximum of each multipole are given in 
Table X and include, as before, the test-mass informa- 
tion. We checked that these fits are compatible with the 
fits given in Table II of Ref. [28]. 

When comparing Fig. 20 with the bottom left panel of 
Fig. 10, we see that, though the effect of having replaced 
by ^A^peak Small, it leads to visible differences. 
In particular, one sees that the frequency evolution near 
merger was more accurately captured in Fig. 10 than in 
Fig. 20. 



Appendix C: Effect of including NQC corrections to 
higher multipoles in the radiation reaction 

In this Appendix we explore the effect of including the 
NQC correction factor in the higher multipole contribu- 
tions to radiation reaction, specifically in some of the 
main subdominant multipoles, /i^/^^, h^3^'~' and h^-^'~' ■ 
[By contrast in the main text we NQC corrected only 
^22^^ in the radiation reaction]. Note that with our 
choice X = oi the argument in pirn{x) we need larger 
NQC modulus correction factors than Ref. [28] which 
used X = ri^/^. Indeed as during the plunge fi^/'^ is 
larger than and as the function pim {x) is a decreasing 
function of its argument, one has, along the EOB dynam- 

ics, (pimiv"^)) > {pim{^'^'^)) ■ Therefore the inclusion 

of NQC corrections for higher multipoles is apriori more 
significant within our EOB setup than within the one of 
Ref. [28] . We focus on the mass ratio g = 6 only, because 
subdominant multipoles do not significantly contribute 
when q ^ 1. Figure 21 compares the phase difference 
for two EOB models: one with the standard /i22-only 
NQC flux correction, and another one wich includes in 
addition the three subleading NQC factors /i^/^'"", h^^^''^ 
and h^^^ ■ One sees that the effect of this inclusion is 
small (actually it is within the NR uncertainties of the 
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g = 6 waveform), but definitely visible. However, the fig- 
ure also shows that by slightly modifying the values of Cg 
(namely using = -122 + 143(1 - 4 x 6/49) = -49.04, 
instead ofag = -122 + 147(1-4x6/49)= -47), one can 
reproduce the previous result. In view of this and in the 
absence of higher accuracy g = 6 NR data, we decided 
to include only the £ = m = 2 NQC correction to the 
radiation reaction. 



Appendix D: Explicit expression of pim{x) and 5tm{y) 

In this Appendix we list the explicit expressions of 
the residual amplitude, pimix), and phase, Sim{y) cor- 
rections that we have implemented in our EOB code. 
They rely on the results of Refs. [5, 51]. We give ex- 
plicit expressions for all multipoles up to £ = 8 included. 
Such expressions are given at the 3"*"^ PN approximation, 
i.e. the 3PN-accurate, ^ results of Ref. [5] are hy- 
bridized with the 5PN-accurate, = 0, terms obtained in 
Ref. [51]. Let us recall that we used here the following 
values of the arguments of these functions: x in 
P£„(a;) and y {H^oB^f^- 
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The "culerlog" functions eulerlog^(x) are defined as 
eulerlog„j(a;) = 7^ + log 2 + i log a; + log m , (D35) 

where is Euler's constant, = 0.577215 . . . and 
log(x) the natural logarithm function. 

Let us now give the explicit expression of the residual 
phase corrections that are implemented in the code. 



For (52m, ^33 and ^31 we list here explicitly both their 
Taylor-expanded (labelled with a "Taylor" superscript) 
and their Fade resummed. The S^m for higher multipoles 
can be given only in Taylor-expanded form and thus the 
label "Taylor" is omitted. The terms in boldface are 
the highest-order known FN terms for v ^ {). They are 
omitted when 7^ 0, and in particular in the computation 
of the Fade approximants, but they arc kept only in the 
computation of the v ~ Q BOB waveform. The Taylor- 
expanded 5hn read 
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(D36) 
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(D40) 
(D41) 
(D42) 
(D43) 
(D44) 
(D45) 



Among these, we used S32, Sim and 6^5 in their Taylor- 
expanded form indicated above. By contrast, for 622, 



I 

S21, S33 and S31 we used (denoting Vy = ,Jy) the following 
Fade resummed expressions (see Section II B 1 for further 
details) 
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